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ABSTRACT 

Two different viscous-inviscid interaction schemes have been 
developed for the analysis of steady, turbulent, transonic, separated 
flows over axisymmetric bodies. The viscous and inviscid solutions are 
coupled through the displacement concept using a transpiration velocity 
approach. In the semi-inverse interaction scheme, the viscous and 
inviscid equations are solved in an explicitly separate manner and the 
displacement thickness distribution is iteratively updated by a simple 
coupling algorithm. In the simultaneous interaction method, local 
solutions of viscous and inviscid equations are treated simultaneously, 
and the displacement thickness is treated as an unknown and is obtained 
as a part of the solution through a global iteration procedure. 

The inviscid flow region is described by a direct finite-difference 
solution of a velocity potential equation in conservative form. The 
potential equation is solved on a numerically generated mesh by an 
^PP^oxinicite factorization (AF2) scheme in the semi-inverse interaction 
method and by a successive line overrelaxation (SI.OR) scheme in the 
simultaneous interaction method. 

The boundary- layer equations are used for the viscous flow region. 
The continuity and momentum equations are solved inversely in a coupled 
manner using a fully implicit finite-difference scheme. The energy 
equation is solved uncoupled. The FLARE approximation is used in the 
reversed flow region and its effectiveness is studied by using a 
windward differencing scheme. 
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The two- layer algebraic turbulence model proposed by Cebeci and 
Smith (1974) and a new one-half equation turbulence model proposed by 
Johnson and King (1985) are utilized to describe the Reynolds stress in 
turbulent flow calculations. Parameters affecting the convergence rate 
of the interaction procedure are discussed. The calculation schemes are 
evaluated by studying 1) an incompressible, laminar, separated flow over 
a flat plate with a trough, 2) a turbulent, transonic, separated flow 
over an axisymmetric boattail with a solid cylindrical plume simulator, 
3) a turbulent, transonic, separated flow over an axisymmetric bump 
attached to a circular cylinder. The predictions are compared with 
experimental data and other available numerical results. The 
simultaneous interaction method becomes more efficient and reliable than 
the semi-inverse method as the separation size grows. The prediction 
obtained by the Johnson-King model is generally in good agreement with 
the measurements, but disagreement is noticeable after the reattachment 
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I . INTRODUCTION 
A. Overview of the Problem 

Encouraging progress has been made in recent years in the 
prediction of complex flow fields. However, much more work needs to be 
done to develop and verify reliable predictive schemes in several areas 
of application. One of the difficult problems in present day 
aerodynamics is the accurate prediction of complex turbulent flows 
occurring in the transonic speed regime. Transonic flows occur in many 
important aerodynamic applications of current technological interest 
including flows over airfoil and engine components. This is largely 
because many modern commercial and military aircraft cruise very 
efficiently at transonic speeds. 

While transonic flow fields contain a variety of interesting and 
unique characteristics, they are especially characterized by the two 
main complicating features of mixed subsonic/supersonic flow and 
substantial viscous effects. In a typical transonic flow, a subsonic 
freestream is accelerated over a convex body surface to form an embedded 
region of supersonic flow adjacent to the body surface. This supersonic 
region is generally terminated by a shock wave that recompresses the 
flow as the flow returns to the subsonic speed. The strength and extent 
of the shock wave increase with freestream Mach number. From a 
mathematical point of view, the transonic flow must be described by a 
nonlinear equation or set of equations of mixed elliptic/hyperbolic 
type, because a subsonic flow region is described by an elliptic 
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equation and supersonic region is described by a hyperbolic equation. 

The boundaries between the elliptic and hyperbolic regions must be 
obtained as a part of the solution. 

Viscous effects are also important in the transonic flow. The 
interaction between a viscous (or boundary) layer and a shock wave is a 
very complicated phenomenon. The principal interaction between the 
shock wave and the boundary layer arises from the displacement thickness 
effect which leads to a thickened effective body causing significant 
changes to the surface pressures and forces. In the inviscid flow 
region, the pressure increases discontinuous ly across the shock wave. 
However, in the inner part of the boundary layer which has subsonic 
velocity, this abrupt increase of the pressure cannot occur. Instead, 
the overall pressure rise takes place over several boundary- layer 
thicknesses. Although the shock wave penetrates the boundary layer and 
generates a significant normal pressure gradient, it is considerably 
weakened and finally vanishes as it reaches the sonic line in the lower 
part of the boundary layer. The flow under the sonic line is retarded 
by this adverse pressure gradient and boundary layer is substantially 
thickened. 

If the strength of retarding influences is sufficiently strong, the 
boundary layer separates from the body surface increasing the streamwise 
spread of the pressure rise. A typical example of the transonic flow 
with the interaction between the normal shock wave and the boundary 
layer is schematically illustrated in Figure 1. When the overall 
pressure rise is sufficient to cause separation, an outgoing weak 
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FIGURE 1. Transonic flow field through the interaction 
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compression wave is generated. Although the compression wave downstream 
of the separation is barely perceptible, it joins the nearly normal 
strong shock some distance outside the boundary layer. To achieve 
continuity of pressure and flow direction downstream of this 
intersection, an incoming weak oblique shock is generated, and a vortex 
sheet with a rapid entropy rise across it runs downstream from the 
intersection. Because of a nonuniform entropy increase across the shock 
wave, the flow outside the boundary layer immediately downstream of the 
weak trailing shock is still supersonic, but with further pressure rise 
as the flow continues downstream the supersonic tongue diminishes and 
finally vanishes (Green, 1970). 

The presence of separation can significantly affect the entire flow 
field. For example, the existence of separation on an airfoil can 
change the lift and drag coefficients considerably so as to degrade the 
control effectiveness of the aircraft. Also, the flow reattachment 
gives rise to heating rates which can far exceed those for an attached 
boundary layer. Therefore, the mechanism of the separation and 
reattachment in transonic flow fields has received a great deal of 
attention in the aerodynamic design process, though it is still far from 
being fully understood. Generally, the process of separation is 
believed to depend only on the properties of the approaching flow 
stream. However, the flow inside the separation bubble is thought to be 
subject to both upstream and downstream influences, thus having elliptic 
characteristics. The size of a separation bubble usually depends on the 
magnitude of the pressure rise, the nature of the disturbance which 
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causes it, the Mach number, and the Reynolds number of the initial flow 
(Green, 1970). However, the size of separation caused by a normal shock 
wave is so influenced by the interaction with the far field that a 
correlation of separation size with the main flow parameters is very 
difficult to obtain. Such an interaction process is complex and our 
present understanding of this phenomenon is very incomplete. 

The complexity of the governing equations that need to be solved to 
predict transonic flows depends on the flow phenomena in question. If 
the shock wave is sufficiently weak that the flow does not separate, the 
viscous effects are often insignificant. In this case, a fairly 
accurate description of the pressure distribution can be obtained from 
solving the equations for inviscid flow only, perhaps even the transonic 
small disturbance (TSD) form of the equations. If the shock wave is 
strong enough to cause separation, the governing equations must include 
viscous effects. The optimum (in terms of cost and accuracy) 
computational strategy for the latter case has not been established. 

The broad classes of competing methods can be labeled as global or 
zonal . 

In the global approach, no assumptions are necessary about the type 
or nature of interaction in the flow field. All regions are computed 
simultaneously with a single set of equations, such as the Navier-Stokes 
equations. The Navier-Stokes equations are generally regarded as the 
basic equations describing most flow phenomena of practical aerodynamic 
interest. The physical phenomena encountered in transonic flows, 
including mixed subsonic/supersonic flow, shock waves, boundary layer, 
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separation and turbulence can be mathematically represented by the time- 
dependent Navier-Stokes equations. The numerical solution of the time- 
dependent Navier-Stokes equations for practical turbulent flow problems 
is, however, presently not feasible with existing computers due to the 
small characteristic length and time scales of the turbulent motion. 

Thus, some form of averaging of the governing equations and turbulence 
modeling are required. Current practice, as well as much research in 
turbulent flows, is based on the use of the time-averaged Navier-Stokes 
(Reynolds) equations with turbulence modeling. Often, the predictions 
based on this approach have failed to predict flow details accurately, 
but shortcomings are frequently attributed to the grid arrangement or 
the turbulence model. 

In the zonal approach, the flow region is divided into subregions 
which have distinct flow characteristics and each subregion is described 
by an appropriately reduced set of governing equations. As in many flow 
situations, the transonic flow field can be divided into the thin shear 
(boundary) layer near solid boundaries and inviscid flow elsewhere, 
assuming that the Reynolds number is sufficiently large. Most commonly, 
the Euler equations or the potential equation have been used for the 
inviscid flow region and either an integral or finite-difference 
representation of the boundary- layer equations for the viscous flow 
region. It is often possible to solve these two sets of equations 
interactively in a consistent manner using what has become known as the 
viscous-inviscid interaction approach. Applications of viscous-inviscid 
interaction schemes based on the above idea to transonic turbulent 
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separated flows can be found in numerous works. Predictions based on 
this approach generally compare well with the solutions based on the 
time-averaged Navier-Stokes equations, but solutions from both 
approaches often exhibit disagreement with experimental measurements for 
transonic separated turbulent flows. These discrepancies in solutions 
based on the viscous-inviscid interaction approach have been attributed 
to inadequate turbulence modeling or errors associated with the standard 
boundary- layer approximation which neglects the normal pressure gradient 
across the boundary layer. 

The accuracy and reliability of turbulent flow predictions are very 
much constrained by the accuracy and generality of the turbulence model 
used to evaluate the Reynolds stresses and heat flux quantities. To a 
large extent, turbulence modeling is the pacing item in the quest for 
improved predictions. Despite considerable research effort, completely 
satisfactory turbulence models have not been identified for many complex 
flows, especially those containing regions of flow reversal. The lack 
of generality is a major shortcoming of turbulence models. Many of the 
studies for viscous transonic flows have used variants of a simple 
algebraic turbulence model (Deiwert, 1976; Baldwin and Lomax, 1978; 
Carter, 1981). Recently, more complex models which solve additional 
differential equations for turbulence parameters have been used. 
Surprisingly, complex models that show better overall predictions in 
many incompressible cases do not provide significant improvement over 
simple algebraic models for transonic flows with large separation. 
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A long range goal of computational fluid dynamics is the 
development of methods which are predictive; that is, methods which can 
be used with reasonable confidence in the absence of confirming 
experimental data. To support this goal, the development of turbulence 
models with reasonable generality could be helpful. Until this has been 
achieved, the identification of models which work well for particular 
classes of flows (with their accompanying range of applicability and 
limitations) would prove useful for design purposes. 

The present study deals with both computational and turbulence 
modeling aspects of predicting transonic flows with strong interaction. 
An objective of the study was to advance the state of the art in 
turbulent flow prediction and to enhance present understanding regarding 
the range of applicability and limitations of computational approaches 
and turbulence models for flows containing separated regions. 

B. Literature Review 

The study of transonic aerodynamics has a long history and a 
considerable amount of research has been done both experimentally and 
theoretically. The present review is intended to provide an indication 
of the state of the art of transonic aerodynamics . It is by no means 
all inclusive. The emphasis was put on works of two-dimensional 
transonic flows with strong interaction between the shock wave and the 
boundary layer, which is closely related to the present study. 
Experimental studies and analytical works are discussed separately for 
convenience. A detailed historical review of transonic flow research 
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can be found in an article by Spreiter (1982). 

1. Experimental work 

Interest in transonic flows was developed as early as the 1900s, 
but realistic experiments only began in the 1940s because of technical 
difficulties in obtaining transonic speeds. However, at present, 
experimental measurements over a variety of flow conditions and body 
configurations are available. Still, only a limited number of 
experimental investigations provide extensive, quantitative measurements 
flow properties in the boundary layer which are necessary for 
understanding the mechanism of the interaction between a shock wave and 
a boundary layer. 

Naturally, most of the transonic flow experiments have focused on 
airfoil-like geometries. Some of the important works are given below. 
Liepmann (1946) performed an experimental study for the transonic flow 
past a 12/6 circular arc airfoil to examine the effect of the boundary 
layer upon the shock wave pattern and pressure distribution. 

Measurements indicated that a change of the boundary layer from the 
laminar to the turbulent regime resulted in a marked change in shock 
wave pattern and surface pressure distribution at the same freestream 
Mach number. Also, the pressure gradient normal to the boundary layer 
was found to be of the same order as the one parallel to the boundary 
layer near the base of the shock wave. 

Measurements of Ackeret et al . (1947) showed that the boundary- 
layer displacement thickness increased rapidly at separation so that the 
surface pressure distribution was significantly modified from what would 
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be predicted from pure inviscid analysis. According to the measurements 
of Seddon (I960), the incident normal shock is bifurcated into a strong 
oblique shock and a weak rear shock due to the viscous interaction. 
Following the trailing weak rear shock, the supersonic stream is 
embedded in the subsonic region just outside the boundary layer. In the 
boundary layer, the flow was shown to undergo successively the three 
processes of shock compression, displacement and rehabilitation. 

Pearcey et al. (1968) carefully examined the scale effects in wind 
tunnel tests. They suggested that the applicability of the flow model 
developed on shock induced separation of turbulent boundary layers on an 
airfoil to full scale behavior could be restricted. Restrictions arise 
because the trailing edge does not include the interaction that 
sometimes occurs between the disturbance at the foot of the shock and a 
subsonic-type rear separation in the continuous adverse gradient further 
downstream . 

McDevitt et al. (1976) performed an experimental study for the 
transonic flow over an 18% thick circular-arc airfoil. By varying the 
peak local Mach number from about 1 to 1.4, both weak and strong 
interactions were observed. Shock induced separation was observed at 
M =0.78, and the effect of changes in Reynolds number on the flow 
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fie^d was appreciable at low Reynolds numbers, but this effect was small 
for Reynolds numbers above 1 x 10^ based on the airfoil chord length. 

For freestream Mach numbers ranging from about 0.76 to 0.78, the airfoil 
flow field was found to be unsteady. Comparisons of the measurements 
and numerical solutions of the Navier -Stokes equations suggested that 


the development of a more accurate turbulence model was necessary when 
the interaction is strong and extensive separation is present. The same 
flow model was tested by Levy (1978) and Seegmiller et al. (1978) at a 
fixed Reynolds number and the same unsteady motion was observed at the 
limited range of freestream Mach numbers indicated above. 

Johnson and Bachalo (1980) reported an experimental study for a 
symmetrical NACA 64A010 airfoil at transonic conditions. Measurements 
were obtained with the freestream Mach number fixed at 0.8 for three 
different angles of attack to vary the intensity of shock wave/boundary 
layer interaction. The effect of varying Reynolds number was found to 
be very small because the transition strip placed at the leading edge 
nullified the effects of natural transition, thus reducing the 
sensitivity to Reynolds number. As the angle of attack was increased, 
the boundary- layer thickness on the airfoil f s upper surface 
significantly increased and the shock wave moved forward. The turbulent 
flow measurements revealed that the turbulence fluctuations attained 
equilibrium with the local mean flow more rapidly than previously 
expected. On this basis, it was suggested that improved turbulence 
modeling was needed at or very near the separation point and that an 
algebraic turbulence model based on the local equilibrium assumption can 
be used downstream of the separation point. 

One of the most heavily investigated transonic flow configurations 
other than the airfoil configuration is the cone-cylinder boattailed 
afterbody. Shrewsbury (1968) studied the effect of boattail juncture 
shape on afterbody drag at transonic speed. He tested eight different 
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afterbody configurations with cylindrical plume simulators for 
freestream Mach numbers ranging from 0.56 to 1.0 and angles of incidence 
ranging from 0 to 8 degrees. The experiment was performed with both 
rounded and sharp sting junctures. A sharp juncture was found to result 
in a slower pressure recovery on the afterbody. This experiment 
included measurements for three-dimensional separated flows. 

Reubush (1974) conducted an experimental investigation for a series 
of eight nacelle-mounted isolated circular-arc boattailed afterbodies 
with both cylindrical plume simulators and real jet exhaust plumes to 
determine the effectiveness of utilizing solid circular cylinders to 
simulate jet exhaust plumes. The experiment was conducted at freestream 
Mach numbers ranging from 0.4 to 1.3 at an angle of attack of 0 degrees 
with various ratios of simulator diameter to nozzle-exit diameter. 
Comparisons of the measurements generally indicated that use of one of 
the larger diameter simulators would approximately result in pressure 
coefficient distributions and drag coefficients of real jet exhaust 
plumes at all Mach numbers. A more detailed discussion of these 
experiments will be given in a later section. 

Later, an investigation was conducted by Reubush and Putnam (1976) 
to determine the effects of variations in Reynolds number on the 
pressure and drag of similar isolated boattails. They found that as the 
Reynolds number was increased, the boattail static pressure coefficients 
in the expansion region of the boattail became more negative, although 
those pressure coefficients in the recompression region of the boattails 
became more positive. These two trends were found to be compensating 
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and there was only a small effect of Reynolds number. Abeyounis (1977) 
measured the separation point using an oil-flow technique for the same 
circular-arc boattails, and suggested a correlation between the 
separation point and the Mach number, but the reattachment point was not 
measured. 

Benek (1979) experimentally investigated transonic flows over two 
different boattail configurations with solid cylindrical plume 
simulators for freestream Mach numbers ranging from 0.6 to 1.3 using a 
laser doppler velocimeter. The two configurations differed in 
smoothness of the boattail so that one produced a fully attached flow 
and the other produced substantial separation. 

There have been several experimental studies conducted on the flow 
field over circular-arc bump attached to a plane wall or cylinder. This 
configuration has the advantage of generating thick, extensive 
separation compared to the airfoil geometry so that more information can 
be obtained about the turbulent boundary- layer flow in the immediate 
vicinity of the separation point. 

Alber et al. (1973) measured the turbulent transonic separated flow 
generated on the rear portion of a two-dimensional circular-arc model 
mounted on the wind tunnel floor using pitot probes. By varying the 
freestream Mach number from 0.3 to 0.8, two different types of 
separation were tested; pressure gradient induced separation and shock 
wave induced separation. At the lower Mach numbers, no shock was 
observed and separation induced by the pressure gradient appeared at the 
aft portion of the bump. As the Mach number was increased up to 1.32, 
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the separation point moved upstream and the reattachment point moved 
downstream. Shock induced separation was observed as the peak local 
Mach number was increased to 1.34. As the Mach number was increased 
further, the shock and the separation and reattachment points moved 
downstream. Velocity profiles downstream of the shock wave were also 
found to be quite different in the pressure gradient induced separation 
and shock induced separation cases. However, the accuracy of the 
measurements of velocity profiles in the reversed flow region is 
uncertain due to possible errors associated with the measuring 
technique . 

Using a Mach Zehnder interferometry, Delery et al. (1976) studied a 
separated flow with strong viscous -inviscid interaction in a two- 
dimensional transonic channel in which a 12% thick circular bump profile 
was placed on the lower channel wall. Their results clearly indicated 
that the interaction with the shock wave brought on a noticeable 
thickening of the boundary layer and a noticeable distortion of the 
velocity distribution. They also found that the appearance of 
separation depended on the shock intensity and the velocity profile 
shape . 

Delery (1983) investigated flows resulting from the interaction of 
shock wave and turbulent boundary layer occurring in a two-dimensional 
transonic channel. The turbulent boundary layer was developed on the 
channel walls and the bumps were placed on both upper and lower channel 
wall to accelerate the flow to supersonic speed. It was found that the 
first part of the interaction process entailed a very large production 
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of turbulence having a very strong anisotropy. In this zone, the normal 
stresses in the momentum and turbulent energy equations were considered 
to be important. He reported that the downstream relaxation toward a 
new equilibrium state was a very gradual process due to the long 
lifetime of the large structures which were formed in the region of 
intense turbulence production. 

Bachalo and Johnson (1979) performed an experimental study for an 
axisymmetric flow model which consisted of an annular circular-arc bump 
affixed to a circular cylinder aligned with the flow direction. 
Measurements were obtained in the NASA Ames 2x2 foot transonic wind 
tunnel at freestream Mach number of 0.875 using a laser velocimeter. 

The separation and reattachment points were determined using oil flow 
visualization. The same bump configuration was tested by Horstman and 
Johnson (1984) at Mach numbers raging from 0.4 to 0.925 in the NASA Ames 
6x6 foot supersonic wind tunnel where the influence of the tunnel 
walls would be much smaller than in the previous smaller wind tunnel. 

For the same test conditions, measurements of the surface pressure 
distribution obtained in these two wind tunnels were almost identical 
except that the shock location obtained in the larger tunnel was always 
about 1% chord length upstream of that in the smaller wind tunnel. As 
the Mach number was increased up to 0.8, the flow remained subcritical, 
but a small region of separation was observed. At Mach number slightly 
over 0.8, a shock wave was formed. As the Mach number was increased 
further, the shock location remained almost the same, but the separated 
region grew rapidly. More detailed observations will be given in a 
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later section. 

2. Analytical work 

a. General Since transonic flow is governed by a nonlinear 

equation or set of equations, analytical studies were very limited in 
the early days of transonic aerodynamics. Therefore, early studies of 
transonic flows depended very heavily on experiments and there were only 
few theoretical results. More recently, due to the general availability 
of digital computers, many numerical studies have been conducted and 
these have accelerated understanding of transonic phenomena. 

Numerical studies of transonic flows can generally be divided into 
two categories depending on whether viscous effects are accounted for or 
not. Early studies were typified by the solution of the inviscid subset 
of the Navier - Stokes equations, such as the transonic small disturbance 
equation solved by Murman and Cole (1971) for two-dimensional flow and 
Bailey and Ballhaus (1972) for three-dimensional flow, the potential 
equation, solved by Jameson (1974), and the Euler equations, solved by 
Magnus and Yoshihara (1970). Comprehensive reviews concerning this 
approach can be found in articles by Yoshihara and Spee (1982), Nixon 
and Kerlick (1982), Holst (1983), and South (1985). These fully 
inviscid analyses often produce approximate descriptions of the 
transonic flow field including the pressure distribution and shock wave 
location. However, they fail to provide an accurate description when 
separation occurs due to shock wave penetration or a strong adverse 
pressure gradient. As observed in many experiments, viscous effects in 
transonic separated flows are so significant that a fully inviscid 


analysis cannot be used to obtain a reasonable solution in such flows. 

The other category consists of numerical methods which include 
viscous effects. Such methods are usually divided according to two 
approaches. The first approach is based on the full set or reduced 
subset of the Navier-Stokes equations and the other is the zonal 
approach based on the coupling of the viscous and inviscid subsets of 
the Navier-Stokes equations (the viscous -inviscid interaction method). 

b. Navier-Stokes solutions-general Examples of the first 
approach which generally uses the time dependent Reynolds-averaged or 
mass-averaged Navier-Stokes equations can be found in numerous studies. 
An informative review on this approach is contained in an article by 
Mehta and Lomax (1982). Some of these studies were devoted to an 
evaluation of the turbulence modeling, generally the algebraic model 
with or without relaxation type modifications or one or two-equation 
models, for shock wave/boundary layer interaction flows. Included in 
such studies are Baldwin et al. (1975), Baldwin and Rose (1975), Deiwert 
(1976), Viegas and Coakley (1977), Levy (1978), Viegas and Horstman 
(1979), and Horstman (1983). Generally, higher-order turbulence models 
such as two-equation models provided a better description of flow 
details than the algebraic turbulence models . However, differences 
between these two models became very small in flows with shock induced 
separation and the overall predictions were in poor agreement with the 
measurements . 

Several studies also utilized a reduced set of equations that falls 


between the full Navier-Stokes equations and the boundary- layer 
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equations. Baldwin and Lomax (1978) applied a thin-layer approximation 
to the Navier-Stokes equations. These equations are somewhat simpler 
than the Navier-Stokes equations, but a substantial amount of computer 
effort is still required to solve them. The parabolized Navier-Stokes 
(PNS) equations have recently gained popularity in many flow 
calculations (Rudman and Rubin, 1968) because of the use of a space 
marching technique, but their application to transonic flow has been 
rare. Recently, however, Khosla and Lai (1984) developed a global PNS 
technique to calculate transonic separated flows. 

c. Navier-Stokes solutions -boattails and bumps A more detailed 

description of previous studies conducted with the above approaches on 
flows of present interest, the boattail afterbody and bump 
configurations, is given below. Holst (1977) and Swanson (1980) solved 
the time-dependent, Reynolds -averaged Navier-Stokes equations using an 
explicit finite-difference method for axisymmetric boattail flows 
measured by Reubush (1974). Both studies used a two-layer algebraic 
turbulence model with a relaxation formula to account for the 
nonequilibrium effects of the interacted turbulent flow. Predictions 
agreed much better with experimental data than results of equilibrium 
turbulence models, but the predictions in the reversed flow region were 
not satisfactory. This suggests that the improvement of the turbulence 
modeling is needed especially in the separated flow region. 

Deiwert (1981) used the thin- layer form of the Reynolds -averaged 
Navier-Stokes equations, developed by Pulliam and Steger (1980), to 
calculate flows over several different axisymmetric boattail 
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configurations experimentally tested by Shrewsbury (1968), Reubush 
(1974) and Benek (1979). A two-layer algebraic turbulence model 
introduced by Baldwin and Lomax (1978) was used. Results agreed well 
qualitatively, but quantitatively some discrepancies were observed in 
the separated flow region; surface pressure was overpredicted in the 
pressure recovery region. These discrepancies were attributed to poor 
turbulence modeling, effects of artificial viscosity, and the influence 
of grid point distribution. 

Khosla and Rubin (1983) applied a composite velocity formulation in 
solving the Navier-Stokes equations for Reubush' s boattail flow. This 
formulation used a multiplicative composite of the appropriate velocity 
representations for the inviscid and viscous flow regions. As a result, 
the equations were structured so that far from the surface of the body 
the continuity equation reduced to the potential equation and the 
momentum equations led to the Bernoulli equation. Swanson et al. (1983) 
extended this method to transonic flow calculations by using the 
artificial compressibility method for embedded supersonic regions. A 
two-layer algebraic turbulence model proposed by Cebeci and Smith (1974) 
was used with the relaxation model suggested by Shang et al. (1976). 
Results showed favorable agreement with experimental data but the 
pressure was still overpredicted in the large separated flow region. 

Solutions of the time-dependent, Reynolds -averaged Navier-Stokes 
equations were obtained by Johnson et al . (1982) for the flow over an 
axisymmetric bump attached to a cylinder, experimentally measured by 
Bachalo and Johnson (1979). The two-equation (k-w) turbulence model 
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proposed by Wilcox and Rubesin (1980) and the algebraic Cebeci-Smith 
turbulence model were used. Even though the two-equation model resulted 
in slightly better predictions than the Cebeci-Smith model, both models 
produced essentially the same results. Particularly, the shock 
locations predicted by both models were in poor agreement with the 
measurements. The authors suggested that this large discrepancy was due 
to the rapid increase in the turbulent shear stresses in the vicinity of 
the shock. 

Using the same Navier-Stokes equations, Horstman and Johnson (1984) 
recalculated this flow for a wider range of freestream Mach numbers and 
compared with their new measurements. At the outer boundary, freestream 
conditions were assumed instead of the inviscid solid wall (slip) 
boundary condition. They also used the two-equation (k-s) model of 
Jones and Launder (1972) with a longitudinal curvature correction. The 
change of the outer boundary condition resulted in a vast improvement in 
calculated results, which indicated that the measurements were free from 
effects of the wind tunnel wall. A slight improvement was observed by 
using the k-£ turbulence model, but the pressure was still overpredicted 
in the recovery region and the separated region was predicted very 
poorly compared to experimental measurements. The computations showed 
no separation for freestream Mach numbers up to 0.8, while the 
experiment showed separation at all Mach numbers tested. For the larger 
freestream Mach numbers, the predictions of the separation points were 
relatively good, while the calculated reattachment points were 
significantly upstream of the measured points. 


21 


Sahu and Danberg (1985) used the thin-layer Navier-Stokes equations 
to predict the transonic separated flow over the same axisymmetric bump 
model. The computations were made with the algebraic turbulence model 
of Baldwin and Lomax (1978) and the two-equation (k-t) turbulence model 
of Chien (1982). The calculated surface pressure distributions with 
both turbulence models were almost identical throughout the interaction 
region. The shock location was well predicted, but the pressure was 
significantly overpredicted downstream of the shock by both models and 
the pressure plateau was not captured by either model. They thought 
that this large discrepancy was due to the poor grid spacing in the 
redeveloping region. But the large error in prediction may also be 
caused by inadequate turbulence modeling. The large difference between 
the predictions of the two turbulence models was noticeable in the 
predicted skin-friction coefficients downstream of the reattachment 
point. The two-equation model resulted in a very strange sharp peak 
downstream of the reattachment point. The predicted reattachment point 
was located about 15% of the bump chord length upstream of the measured 
point. Therefore, the accuracy of the k-E turbulence model of Chien 
(1982) is somewhat questionable for this flow. 

The Johnson-King turbulence model was used by Johnson (1985) in 
solutions to the Navier-Stokes equations for axisymmetric bump flows of 
Bachalo and Johnson (1979). Calculated results were compared with the 
measurements for a wide range of freestream Mach numbers and were 
generally observed to be in very good agreement with the measurements. 
However, a slight overprediction of the displacement thickness was 
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noticeable downstream of bump-sting juncture, and a slow recovery of the 
flow toward the equilibrium state was evident downstream of the 
reattachment point. Still, the Johnson-King turbulence model was found 
to provide very good overall predictions for flows with strong 
interaction. The model requires very little more computational effort 
than the equilibrium, algebraic models. 

d. Viscous-inviscid interaction-general Poor predictions by 

the fully inviscid numerical schemes and the high computing cost of 
Navier-Stokes solutions have accelerated the development of the viscous- 
inviscid interaction approach in recent years. A more comprehensive 
discussion of this approach can be found in the recent reviews by 
Le Balleur (1981a), Lock (1981), and Melnik (1981). This approach 
generally requires the iterative matching of the viscous and inviscid 
solutions through displacement thickness coupling. In the initial 
development stage, the classical direct matching method was used as can 
be found in the works of Bauer et al. (1975), Bavitz (1975), and Collyer 
and Lock (1979). The direct interaction method cannot be applied to 
separated flows because the boundary- layer solution breaks down at the 
point of separation. In order to overcome this singularity problem, 
several asymptotic and empirical ideas were developed (Reubush and 
Putnam, 1976; Wilmoth, 1977). However, these methods could not provide 
significantly accurate solutions due to built-in empiricism. 

The separation singularity problem in the boundary- layer equations 
can be eliminated by the use of an inverse method. With the inverse 
method, prescribing the displacement thickness (Catherall and Mangier, 



1966) or skin friction (Klineberg and Steger, 1974), the boundary- layer 
solution becomes regular. However, the displacement thickness or skin- 
friction distributions are not known and must be obtained as a part of 
solution if the method is to be predictive. Carter and Wornom (1975) 
coupled an inverse boundary- layer solution with an inverse solution of 
Cauchy integral equation for the incompressible flow calculations. This 
approach has received little support because of its slow convergence 
rate. Also, it would seem difficult to implement an inverse formulation 
for transonic inviscid flow. 

Le Balleur (1978) introduced a semi-inverse interaction method for 
the calculation of transonic separated flow in which an inverse solution 
to the integral boundary- layer equations was coupled to a direct 
inviscid solution. A similar approach was developed by Carter (1979) 
and Kwon and Pletcher (1979) in their incompressible separated flow 
analyses, except that the boundary- layer equations were solved in a 
finite-difference form. This approach was also extended to calculations 
of transonic flows in many applications. 

The integral form of the boundary- layer equations has been used 
with various inviscid flow formulations. Whitfield et al. (1981) and 
Murman and Bussing (1984) used the Euler equations; Lee and Van Dalsem 
(1981) used the full potential equation; Melnik and Brook (1985) used 
the full potential equation with an entropy correction. The interaction 
schemes of Carter (1981), Van Dalsem and Steger (1983), and Murphy and 
King (1983) used the finite-difference form of the boundary- layer 
equations for the viscous region and the full potential equation for the 
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inviscid region. Although these schemes gave reasonable predictions for 
some transonic flows, discrepancies between calculated and measured 
surface pressure distributions were often observed in separated flow 
regions . 

In a series of numerical studies by Carter and Hafez (1982) and 
Carter et al. (1983), attempts were made to account for normal pressure 
gradients and the effects of embedded shocks. A compressible stream 
function formulation was used for the inviscid flow to take into account 
the rotational flow effects in the outer region of the boundary layer 
downstream of the shock. The results indicated that for transonic shock 
induced separation, the effects of displacement thickness interaction 
dominated over those produced by embedded shocks and normal pressure 
gradients. It was concluded that the correct turbulence modeling is 
more important in obtaining good predictions than inclusion of normal 
pressure gradients. 

In another approach to the interaction problems known as the 
simultaneous interaction method, local solutions for the pressure and 
displacement thickness are obtained simultaneously so that they can 
mutually satisfy the local viscous and inviscid relations. Most of the 
early applications of this approach were in fully supersonic or subsonic 
flow cases where the outer inviscid flow can be described by a simple 
linear theory or an integral relation (Lees and Reeves, 1964; Crimi and 
Reeves, 1976; Veldman, 1979, 1981; Davis and Werle, 1981; Davis, 1984). 
Moses et al. (1978), on the other hand, used an integral method for the 
boundary- layer equations and a finite-difference method for the solution 
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of the Laplace equation for stream function in the inviscid flow region 
for a laminar incompressible flow. 

A similar approach using the integral method for the boundary- layer 
equations was then extended to steady transonic flows by Wai and 
Yoshihara (1981) and to unsteady problems by Houwink and Veldman (1984) 
and Le Balleur (1984). A fully finite-difference simultaneous 
interaction method was first introduced by Edwards and Carter (1985) for 
laminar incompressible separated calculations. In some of the above 
studies based on the simultaneous interaction approach, the simultaneous 
solution procedure is followed by the pure inviscid or viscous 
calculations (Veldman, 1981; Veldman and Lindhout, 1983; Edwards and 
Carter, 1985). This approach is often called a quasi-simultaneous 
interaction method. 

e. Viscous-inviscid interact ion-boattails and bumps The 

following is a more detailed review of calculation results obtained 
using the viscous-inviscid interaction approach for the boattail and 
bump flows. Chow et al. (1975) applied a conventional direct viscous- 
inviscid interaction approach to calculate the transonic flow over 
axisymmetric boattail bodies tested by Reubush (1974). The full 
potential equation was solved for the inviscid flow region using the 
line relaxation scheme of South and Jameson (1973). The direct solution 
of the integral boundary- layer equations was used for the viscous flow. 
Inviscid solutions were obtained over the equivalent body which was 
obtained by adding the displacement thickness to the original body to 
account for the viscous effects. The displacement thickness was 
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adjusted using a underrelaxation factor of 2/3 after each solution sweep 
of the viscous and inviscid equations. The full potential equation was 
found to yield better agreement with the experimental results than the 
transonic small disturbance equation. However, the direct calculation 
method was not applicable to separated flow. 

Reubush and Putnam (1976) also employed a direct interaction method 
to calculate the separated flow over the same axisymmetric boattail 
configuration. The inviscid flow region was calculated by a linearized 
potential flow (panel) method, developed by Hess and Smith (1967) for 
incompressible flows, with a compressibility correction. Therefore, 
this method was limited to a fully subsonic flow. The viscous flow 
region was described by the integral formulation of the boundary- layer 
equations developed by Reshotko and Tucker (1957). To avoid the 
singularity problem occurring at the point of separation, a 
discriminating streamline method, developed by Presz (1974), was used. 
This discriminating streamline, which separates the reversed flow region 
from the outer boundary layer so that the boundary layer is treated as 
fully attached, was also used as an effective solid body surface. The 
discriminating streamline was given as a conical surface diverging from 
the model surface at an angle which was dependent primarily on the local 
Mach number at the point of separation. A straight line curve fit to 
the experimental data of Presz was used to determine the divergence 
angle. The separation point, which was extremely critical to the 
solution, was determined from the empirical relation of Page (1961) . 

Then the viscous and inviscid flow regions were solved sequentially in a 
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direct mode to update the displacement thickness and surface pressure 
distribution until convergence was obtained. Results showed reasonable 
agreement with experimental data, but still large discrepancies were 
evident in flows with higher Mach numbers. This method is not regarded 
as a truly predictive one because the divergence angle is entirely 
dependent on the experimental data and the separation point is 
determined empirically. A similar approach was used by Presz et al. 
(1978). 

The above method was extended by Wilmoth (1977) to transonic 
calculations by using the line relaxation method of South and Jameson 
(1973) for the finite-difference solution of the full potential 
equation. Also, the discriminating streamline was given as a straight 
line connecting the predetermined separation and reattachment points. 
The separation point was taken from the oil flow measurements of 
Abeyounis (1977) and the reattachment point was assumed to be at the 
point of maximum surface pressure. Therefore, this method cannot be 
regarded as a truly predictive one, either. Better predictions were 
obtained with the experimentally determined discriminating streamline 
than with the empirically determined one. Also, it was found that the 
prediction was very sensitive to the shape and location of 
discriminating streamline. This suggests the need for a truly 
predictive method in which the separated region can be found as a part 
of the solution. Results also indicated that the conservative 
differencing of the potential equation might be necessary when viscous 
effects were included. A similar direct interaction approach was used 
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by Cosner and Bower (1977), but predictions in flows with significant 
separation were poor. 

Dash and Pergament (1978) used a finite-difference formulation of 

the boundary- layer equations in their direct interaction scheme. The 

inviscid region was also described by the full potential equation. 

Instead of a solid plume simulator, a jet entrainment correction was 

used to yield an effective plume boundary. They employed the two- layer 
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algebraic and two-equation (k-E ) turbulence models. The k-E model 
provided better overall results. However, no consideration was given to 
solving separated flow. 

Kuhn (1980) applied the semi-inverse viscous-inviscid interaction 
approach to calculate the same boattail flows. The relaxation method 
used by Wilmoth (1977) was used for the direct solution of the inviscid 
flow equations. The integral form of the boundary- layer equations was 
solved in an inverse mode for the reversed flow region. The 
displacement surface in the reversed flow region, which was prescribed 
as a boundary condition for the boundary- layer equations and as an 
equivalent body surface for the inviscid calculation, was assumed to be 
conical starting at the separation point. The inviscid and viscous 
solutions were alternatively repeated until convergence was obtained. 
The divergence angle and separation point were adjusted in a way to 
reduce the root mean square error between the inviscid and viscous edge 
velocities in the iteration process. The calculated results generally 
agreed well with experimental data. Due to the assumption of linear 
displacement thickness profile, this method cannot generally be applied 
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accurately to other cases such as the separated flow over a flat plate. 

A semi-inverse interaction method was also used by Carter and Vasta 
(1982a) to calculate the axisymmetric boattail flows experimentally 
studied by Reubush (1974) and the axisymmetric bump flow measured by 
Bachalo and Johnson (1979) . The conservative full potential equation 
was solved by a line relaxation scheme in the inviscid flow region. 
Unlike the method of Kuhn (1980), the finite-difference method was used 
to obtain the inverse solution of the boundary- layer equations. The 
standard first-order boundary- layer analysis was also modified to 
account for the effects of the normal pressure gradient which may be 
significant in strongly interacting flows. The viscous and inviscid 
calculations were coupled using a transpiration velocity boundary 
condition for the inviscid flow that was related to the streamwise 
gradient of the displacement thickness following Lighthill (1958). The 
displacement thickness was updated with a formula suggested by Carter 
(1979), which was based on the mismatch of the inviscid and viscous edge 
velocities. With an algebraic turbulence model, the pressure was 
overpredicted and the separation point was too far downstream. Using 
the relaxation model of Shang and Hankey (1975), better predictions were 
observed. Still, predictions in the reversed region appeared to need 
improvement. The effects of normal pressure gradient were not believed 
to be of importance for this flow. 

Whitfield et al . (1981) calculated the transonic interacting flow 
over a planar bump, experimentally studied by Altstatt (1977), using the 
semi-inverse viscous-inviscid interaction method. The time-dependent 
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Euler equations in conservative form were solved in a finite-volume 
formulation for the inviscid flow calculation. The boundary- layer 
equations were solved by an inverse integral technique. The viscous and 
inviscid calculations were coupled using a transpiration velocity. The 
update formula for displacement thickness was identical to that used by 
Carter (1979). Reasonably good agreement between calculated and 
measured surface pressure data was observed. No specific comparisons 
for other quantities were available. 

A semi-inverse interaction method was used by Carter (1981) to 
calculate the transonic flow over the axisymmetric bump experimentally 
studied by Bachalo and Johnson (1979). Carter's scheme utilized a 
finite-difference inverse procedure for the boundary- layer equations. 

The full potential equation was solved by a line relaxation technique. 
The same coupling procedure proposed by Carter (1979) for incompressible 
flow was used. Predictions obtained with the two- layer algebraic 
Cebeci-Smith model showed poor agreement with the measurements; the 
pressure was significantly overpredicted in the vicinity of the bump- 
sting juncture and the separation point was predicted too far 
downstream. The discrepancies were attributed to inadequate turbulence 
modeling. With a reduced value of the Clauser constant in the Cebeci- 
Smith turbulence model, better predictions were observed. The above 
studies based on the viscous -inviscid interaction method except for the 


direct method are summarized in Table 1. 


31 


TABLE 1. The summary of previous studies based on the viscous-inviscid 
interaction method 


References 

Inter- 

action 

Viscous 

analysis 

Inviscid 

analysis 

Appli- 

cation 

Carter and Wornom (1975) 

IV 

FD 

IN 

SB 

Le Balleur (1978) 

SE 

IN 

FD 

TR 

Hoses et al. (1978) 

SI 

IN 

FD 

SB 

Carter (1979) 

SE 

FD 

FD 

SB 

Kwon and Pletcher (1979) 

SE 

FD 

IN 

SB 

Veldman (1979) 

QS 

FD 

IN 

SB 

Kuhn (1980) 

SE 

FD 

FD 

TR 

Whitfield et al. (1981) 

SE 

IN 

FD 

TR 

Lee and Van Dalsem (1981) 

SE 

IN 

FD 

TR 

Carter (1981) 

SE 

FD 

FD 

TR 

Kwon and Pletcher (1981) 

SE 

FD 

FD 

SB 

Veldman (1981) 

QS 

FD 

IN 

SB 

Davis and Verle (1981) 

QS 

FD 

IN 

SB 

Wai and Yoshihara (1981) 

SI 

IN 

FD 

TR 

Carter and Hafez (1982) 

SE 

FD 

FD 

TR 

Carter and Vasta (1982b) 

SE 

FD 

IN 

SB 

Van Dalsem and Steger (1983) 

SE 

FD 

FD 

TR 

Murphy and King (1983) 

SE 

FD 

FD 

TR 

Carter et al. (1983) 

SE 

FD 

FD 

TR 

Veldman and Lindhout (1983) 

QS 

FD 

IN 

SB 

Davis (1984) 

SI 

FD 

IN 

SB , SP 

Houwink and Veldman (1984) 

SI 

IN 

FD 

TR 

Le Balleur (1984) 

SI 

IN 

FD 

TR 

Melnik and Brook (1985) 

SE 

IN 

FD 

TR 

Edwards and Carter (1985) 

QS 

FD 

FD 

SB 


IV 

inverse 

IN : 

integral 

SB 

subsonic 

SE 

semi-inverse 

FD : 

f inite-dif ference 

SP 

supersonic 

QS 

quasi-simultaneous 



TR 

transonic 

SI 

simultaneous 
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C. Scope and Contributions of the Present Study 

Previous studies of transonic flows with strong viscous -inviscid 
interaction described above have used either the time-averaged Navier- 
Stokes equations or the zonal (viscous-inviscid interaction) approach. 
Generally, the Navier-Stokes solutions provide a reasonably good 
predictions for flows with strong interaction between a shock wave and a 
boundary layer, and separated flows can be handled without special 
treatments. However, the numerical solution of Navier-Stokes equations 
still requires a large computational effort, especially for turbulent 
flows. Fine grids are often needed to obtain reasonable solutions 
because the use of wall functions for compressible turbulent flows is 
not well established. Therefore, turbulent flow calculations usually 
require large computing times or they exhibit relatively large errors 
caused by the use of a coarse mesh in order to reduce the computational 
effort. This situation is expected to improve with time, however, due 
to continuing improvements in algorithms and computing machines. 

The need to improve computational efficiency for transonic flows 
provided motivation for the present study. Viscous-inviscid interaction 
methods have been used successfully in the calculation of incompressible 
separated flows by many investigators. As discussed in the preceding 
section, predictions obtained by the interaction method showed 
reasonably good agreement with solutions of the Navier-Stokes equations 
for strongly interacting transonic flows. 

The objective of the present study was to develop an efficient and 
robust viscous-inviscid interaction scheme that can be used to predict 


transonic flows with strong interaction between a shock wave and a 
boundary layer. First, the limit of applicability and efficiency of the 
semi- inverse interaction scheme, developed by Kwon and Pletcher (1979) 
and Carter (1979) was studied. Then a new interaction method was sought 
to enhance the efficiency and applicability of the overall interaction 
approach . 

In the present semi-inverse interaction method, the flow domain was 
divided into the irrotational inviscid and viscous regions. The 
inviscid flow region was described by a direct finite-difference 
solution of a full potential equation in conservative form. The 
inviscid solution procedure was based on the iterative approximate 
factorization (AF2) algorithm which was implemented into the computer 
code, TAIR, by Dougherty et al. (1981). The TAIR code was modified in a 
manner that permitted calculation of any specified portion of an 
axisymmetric flow field having nonperiodic boundary conditions. The 
boundary- layer equations were solved by an inverse finite-difference 
method in the viscous flow region. In order to increase the efficiency 
of the overall interaction scheme and eliminate the possible stagnation 
point problem of the inviscid solution, a shear-layer coordinate system 
was adopted as suggested by Carter (1981). Its influence on the final 
solution and the global convergence behavior was evaluated by using 
several different shear- layer coordinates. 

The viscous solution was coupled with the inviscid solution using 
the transpiration velocity formulation to account for the displacement 
effect. The updating algorithm of Carter (1981), which was used to 
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update the displacement thickness distribution through an iterative 
procedure based on the mismatch between inviscid and viscous edge 
velocities, was modified to increase the convergence rate of the global 
iteration procedure. However, this type of coupling algorithm has been 
known to be unstable for extensively separated flow regions in a 
supersonic stream as suggested by Wigton and Holt (1981). This trend 
has been experienced in the present study. Therefore, the development 
of a more robust interaction scheme was needed for the prediction of 
general separated transonic flows with strong interaction. 

A new simultaneous viscous-inviscid interaction method was 
developed in the present study. In this new interaction coupling 
algorithm, the inviscid and viscous solutions were obtained 
simultaneously in order to eliminate numerical problems associated with 
the semi-inverse interaction method in the calculation of extensively 
separated flow. The simultaneous solution procedure was made possible 
through a localized implicit treatment of two sets of equations. 
Therefore, the updated displacement thickness required for a better 
agreement between the inviscid and viscous pressure distributions was 
obtained as a part of the iterative solution procedure, not with a 
separate interaction algorithm as used in the semi-inverse interaction 
method. To simplify the local treatment of the inviscid solution 
procedure, the conservative full potential equation was solved by a 
simple successive line overrelaxation (SLOR) scheme. 

The solution procedure for the boundary- layer equations was almost 
identical to that of the semi-inverse method except that the 


displacement thickness was treated as unknown to be obtained as a part 
of solution. It was found that most of the computing time was taken by 
the Newton linearization procedure in the boundary- layer equations. A 
pseudo-time dependent approach was developed to reduce the computing 
time needed in the linearization procedure. Since the separation was 
very large in several test cases, there was some concern about the 
marching procedure of the finite-difference scheme in the reversed flow 
region. To verify the effectiveness of the FLARE approximation, a 
windward differencing scheme was also employed in the separated flow 
region for some of the calculations for the purpose of comparison. The 
multiple sweep procedure required for the windward differencing was not 
used separately but was achieved by the global interaction iterative 
procedure. To the authors 1 knowledge, the present new interaction 
scheme is the first simultaneous interaction scheme to be applied to the 
transonic flow regime employing the finite-difference method for both 
the inviscid and viscous equations. 

These two viscous - inviscid interaction methods were demonstrated by 
computations of several different flows including the two-dimensional 
laminar incompressible separated flow over a flat plate with a trough, 
first numerically studied by Carter and Wornom (1975), turbulent 
transonic flow over an axisymmetric boattail afterbody with a solid 
cylindrical plume simulator experimentally tested by Reubush (1974), and 
the turbulent transonic flow over an axisymmetric bump attached to a 
cylinder experimentally measured by Bachalo and Johnson (1979), 
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In the calculation of turbulent transonic flow, turbulence modeling 
is of great importance if separation exists. In the previous studies 
employing either the viscous-inviscid interaction method or the Navier- 
Stokes equations, predictions obtained by both algebraic and two- 
equation turbulence models were generally in poor agreement with 
experimental data, specially in the shock induced separated region. 

These discrepancies are believed to be largely due to shortcomings of 
turbulence modeling. 

The calculated results also suggested that the nonequilibrium 
effects must be accounted for correctly. Johnson and King (1985) 
recently proposed a new turbulence model designed especially for 
turbulent boundary layers in strong adverse pressure gradient with 
separation. This model makes use of an ordinary differential equation 
for the maximum Reynolds shear stress which provides a velocity scale 
for turbulent viscosity. In the present study, this new Johnson-King 
turbulence model was further evaluated in comparison with the algebraic 
Cebeci-Smith model. Observations are made about the effect of varying 
parameters in the Johnson-King model and the influence of the choice of 
turbulence model on the convergence properties of the numerical scheme. 

The contributions of the present study can be summarized as 
follows : 

1. A modification to the semi-inverse interaction method was developed 
to improve the convergence speed of the global iteration procedure. 

2. A shear-layer coordinate system was adopted and its effects upon 
final solutions and convergence were evaluated. 
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3. The validity of the FLARE approximation was evaluated in the shock 
induced separation region by comparing with results obtained with the 
windward differencing. 

4. The pseudo-time dependent approach was developed to reduce the 
computing time needed for the Newton linearization procedure when the 
multiple sweeps of the boundary- layer equations were used. 

5. A new simultaneous interaction method based on the fully finite- 
difference scheme was developed for transonic flow to provide a more 
robust and efficient viscous -inviscid interaction algorithm than the 
semi-inverse interaction scheme. 

6. The Johnson-King model was evaluated in a fully predictive viscous - 
inviscid interaction calculation scheme. 

7. The effect of the choice of turbulence model on the interaction 
convergence behavior was observed. 

In the following chapter, the basic conservation laws needed to 
derive the governing equations for transonic flows are presented. Also 
in Chapter II, the viscous analysis based on the boundary- layer 
equations is presented together with the discussion of turbulence 
modeling. Next, the inviscid analysis based on the full potential 
equation formulation is discussed in Chapter III. In Chapter IV, the 
present viscous - inviscid interaction methods are discussed. Finally, 
Chapter V presents the calculation results in comparison with the 
available experimental measurements and other numerical solutions, which 
is followed by the conclusions in Chapter VI. 
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II. VISCOUS ANALYSIS 

In this chapter, the detailed description of the solution procedure 
for viscous flows will be presented. The general conservation 
principles are described first. The governing equations for 
compressible turbulent boundary- layer flows are then developed from 
those general conservation statements. Turbulence models for both 
equilibrium and nonequilibrium flows are discussed. The numerical 
methods used to solve the equations are also described. 

A. Laws of Conservation 

The fundamental equations of fluid dynamics are based on three 
universal laws of conservation of mass, momentum and energy. The 
governing equations are derived by applying these conservation laws to a 
uniform, homogeneous fluid without mass-diffusion or finite-rate 
chemical reactions. 

1. Conservation of mass 

Using the Eulerian approach, conservation of mass applied to a 
fluid passing through an infinitesimal, fixed control volume yields the 
continuity equation which can be written in vector notation as 

+ V *(PV) = 0 (2.1) 

where p is the fluid density and V is the fluid velocity vector. 
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2. Conservation of momentum 

Conservation of momentum, which is an application of Newton's 
second law to a fixed control volume, can be expressed as 


f-(pV) + V»pVV = pf + V*ff 

0 1 


( 2 . 2 ) 


— > ■ — ^ 

where f is a body force and II is a stress tensor. With a constitutive 

relation for an isotropic, Newtonian fluid based on Stoke' s hypothesis, 

— ► 

the stress tensor, II, is reduced to 


R.. = - p6.. + t.. 

ij ij 


where 6 . is the Kronecker delta function ( 6 . . = 1 if i 
ij 1 J 

if i = i ) and t . . is a viscous stress tensor given by 

j ij 


(2.3) 


i and 6 . . = 0 
ij 


3u . 3u . 3u 

’ij - » + ’ I 5 ij 5*1 


(2.4) 


for Cartesian coordinates. 

Substituting Equation (2.3) into Equation (2.2) and using Equation 
(2.1), the well known Navier-Stokes equation is obtained as 


3V -» -* _ „ -*• 

pf^ + pV-VV = pf - Vp + V* t 

o t 


(2.5) 


3. Conservation of energy 

The conservation law for energy is a statement of the First Law of 
Thermodynamics for fluid passing through a fixed control volume. This 
yields the energy equation written in terms of total enthalpy, H, as 
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pff + P V*VH = |® + V(x*V - q) 
H is defined as 


( 2 . 6 ) 


H = h + | V*V (2.7) 

where h is enthalpy. Using Fourier's law for heat transfer by 
conduction, the heat flux vector, q, can be expressed as 

q = - kVT ( 2 . 8 ) 

where k is the thermal conductivity. 

4. Equation of state 

In order to close the above set of conservation equations, the 
relationships between thermodynamic properties and the transport 
properties must be specified. Since the fluid of interest is air at 
transonic speed, the perfect gas equation of state can be applied with 
little error 


p - pRT (2.9) 

where R is the gas constant. Also for a perfect gas, the following 
relationships hold: 


h = c T 
P 


C 

% 


R 

c = v - T 

v 3f- 1 


*R 

C p Y-l 


( 2 . 10 ) 


where c^ is the specific heat at constant volume, c^ is the specific 
heat at constant pressure and % is the ratio of specific heats. 
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For the coefficient of viscosity, Sutherland's equation was used, 


3/2 


V = c 


1T + C, 


where C and are constants for a given gas. The thermal 
conductivity, k, was evaluated by 


( 2 . 11 ) 


C V 

k = — jj- 
Pr 


( 2 . 12 ) 


where Pr is the Prandtl number, which was assumed to be constant. 


B. Turbulent Flows 


Most flows occurring in nature and in practical applications are 
turbulent. The scientific study of turbulent flow spans approximately 
one century and has resulted in significant progress in many directions. 
Our understanding of turbulent flow is, however, very incomplete. 

Turbulent fluid motion is defined by Hinze (1975) as an irregular 
condition of flow in which the various quantities show a random 
variation in time and space coordinates, so that statistically distinct 
average values can be discerned, and is often characterized with a wide 
range of frequencies and length scales. The size of the largest scale 
is determined mainly by the characteristic dimension of the main flow, 
while the size of the smallest is determined by the fluid viscosity. 

The large scale motion is believed to carry most of the energy and 
momentum in turbulence. The energy is continuously transferred from the 
largest through the intermediate to the smallest scales, where the 
energy is dissipated as heat (Reynolds, 1974). 


i 


/ « 
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The unsteady full Navier-Stokes equations are generally considered 
to describe turbulent flows as well as laminar flows. When the Navier- 
Stokes equations are used to obtain the solution for fluid motion, 
numerical methods must be used instead of analytical procedures because 
of the highly nonlinear characteristics of the equations. However, time 
and space scales of turbulent motion are so small that the large number 
°f grid points and the small size of the time steps required to 
discretize the equations for computer simulation are beyond the 
capabilities of present day computers at least for practical problems. 

The most common practice in the computation of turbulent flows at 
the present time is to solve the time-averaged Navier-Stokes equations, 
which are often referred to as the Reynolds equations of motion, in 
place of the instantaneous equations. These time-averaged equations are 
derived by replacing the instantaneous quantities by the sum of their 
time-mean and fluctuating quantities as 


f = f + f" 


(2.13) 


Time-averaged quantities are denoted by overbars and are defined 


as 


_ ^ t+At 

f = H 1 f(t) dt 


(2.14) 


where At is large compared to the time scale of random fluctuations 
associated with the turbulence but smaller than the time scale of 
unsteady mean motion. For a fluctuating quantity, the time average is 
zero by definition. 


For the treatment of compressible flows and mixtures of gases, it 
is convenient to use the mass -weighted averaging suggested by Favre 
(1965), which removes density fluctuations from the time-averaged 
equations of motion. A mass -weighted average is defined by 

f = f + f' (2.15) 


where 



P 


Hereafter, steady flow is assumed and body forces are neglected. 
Using Equation (2.15) for the velocity and enthalpy and Equation (2.13) 
for the density and pressure, the instantaneous quantities in Equations 
(2.1), (2.5), and (2.6) are replaced with their mean and fluctuating 
quantities. Using the time-averaging procedure and canceling vanishing 
terms in those equations, the mean conservation equations for mass, 
momentum, and energy are obtained as 
continuity 


V*(pV) = 0 


(2.17) 


momentum 


energy 


V(pVV) = - Vp + V(T - pW) 


(2.18) 
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C. Coordinate System 

The coordinate system chosen for the present analysis is presented 
in Figure 2. In the case of axisymmetric external flows, it is common 
to use axisymmetric body intrinsic curvilinear coordinates where the 
abscissa, x, is measured along the body surface and the ordinate, y, is 
measured normal to the body surface as shown in Figure 2. Also, the 
radial distance from the axis of symmetry is noted by r. 

The boundary- layer approximation, which will be discussed in detail 
in the following section, is based on the assumption that the dominant 
viscous shear force is parallel to the body surface. However, when the 
flow separates, the flow direction may deviate substantially from the 
direction determined by the tangent to the body surface. If the 
separated flow region is large, the validity of the boundary- layer 
approximation might be questioned. 

Werle and Verdon (1980) proposed the use of a shear-layer 
coordinate system, shown in Figure 3, for viscous separated flows over a 
blunt trailing edge. This coordinate system is chosen so as to align 
with the predominant direction of the separated viscous flow by adding 
the thickness t to the original body coordinate as shown in Figure 3. 

The boundary- layer approximation is then applied with respect to this 
coordinate system in place of a body-oriented coordinate system. 

Another important advantage of using this shear-layer coordinate system 
arises in the calculation of the inviscid flow in a viscous-inviscid 
interaction procedure. For a body shape with a slope discontinuity, the 
inviscid flow is solved over the shear- layer coordinate by adding the 
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FIGURE 2. Coordinate system for the viscous analysis 
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thickness t to the actual body surface, which eliminates artificial 
inviscid stagnation points thus simplifying the inviscid calculation as 
well as the overall interaction procedure. 

As pointed out by Carter (1981), corrections must be made in the 
viscous-inviscid interaction procedure to account for the change in 
coordinates. The choice of the shear- layer coordinate direction is 
arbitrary but it should approximately align with the main direction of 
flow in viscous regions. Further details on the shear- layer coordinate 
system will be presented in a later section. Calculations were also 
made using the body-oriented coordinate system indicated by Figure 2 in 
the present study. 


D. Boundary-Layer Approximation 

It is well known from experimental observations that, at large 
Reynolds numbers, the effects of viscosity become increasingly confined 
to a narrow region near a solid boundary. In such flow regions, the 
governing equations can be simplified considerably by using Prandtl’s 
boundary- layer approximation. Assuming that the viscous layer is thin 
relative to the characteristic dimension of the object immersed in the 
flow and that the largest viscous term is of the same approximate 
magnitude as any inertia term, an order of magnitude analysis can be 
used to obtain a simplified set of equations, i.e., boundary- layer 
equations . 

This order of magnitude reduction for compressible turbulent flow 
is so lengthy that only the important points will be discussed here (see 


Cebeci and Smith (1974)). The analysis is based on the assumption: 


3 _ 

3y 



u » v 


f » f' 


( 2 . 20 ) 


To estimate the magnitude of turbulence quantities, experimental 
observations must be used. For compressible flows, the magnitude of 
fluctuating terms of density and pressure in addition to velocity and 
temperature must be estimated. If fluctuations of the pressure and 
total temperature, T , are negligibly small, i.e., 

, T 


-£- « 1 Tjr « 1 

P T 

o 

where T q = T + 0.5 u 2 /c^, then temperature fluctuations 
isobaric (Bradshaw, 1977) as 


( 2 . 21 ) 


are nearly 


^ - IL - (x n «2 u ’ 

P T ~ S' 


( 2 . 22 ) 


In the region near the wall, velocity fluctuations are relatively large, 
but the Mach number M is usually small. In the outer region where the 
Mach number is large, velocity fluctuations are generally small. 
Therefore, density fluctuations are small everywhere across the boundary 
layer for Mach numbers even up to 5. Even though there is evidence 
(Kistler and Chen, 1963) that the magnitude of the pressure fluctuation 
is appreciable at Mach number around 5 and is expected to increase with 
increasing Mach number, the pressure fluctuations are usually assumed to 
be negligible and this seems to be a very sound assumption for transonic 
flows. The effects of compressibility upon fluctuations of transport 
properties like the viscosity, conductivity and specific heats are not 
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well understood, but those fluctuations can usually be neglected 
(Bradshaw and Ferriss, 1971). Consequently, the turbulence structure 
of the boundary layer for Mach number up through 5 is generally believed 
to be the same as in the low speed flow, that is, incompressible flow. 
Also, under the boundary- layer approximation, the following relations 
between the time-mean and mass-weighted mean values are found (Cebeci 
and Smith, 1974) 

u = u "h = H (2.23) 


1. Governing equations 

With the above boundary- layer approximation, Equations 
(2 . 17) -(2 . 19) are reduced to the following boundary- layer equations for 
axisymmetric compressible turbulent flows: 
continuity 


3 3 — 

-(pur) + ^(pvr) = 0 

momentum 


(2.24) 


— 3u dp 1 3 r , — 3u — i 

pu 3^ + p v ^ = ' ^ + 7 - pu v )] 


3y 


dx 


r 3y 1 3y 


(2.25) 


energy 


— 3H , — 3H 1 3 . ,p 3H 
pu ta + ■ T a5 (r( p; 


- pH' v' + 


,, 1 . — 3u. 

(1 ‘ p^V’ 


(2.26) 


where r is set to unity for a plane two-dimensional case. These 
equations are also valid for laminar flows when the terms involving the 
fluctuating quantities are set to zero. 
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In obtaining the boundary- layer equations, second derivatives of 
the velocity component in the streamwise, x, direction have been 
neglected along with the entire momentum equation for the transverse, y, 
direction. As a result, the pressure gradient in the transverse 
direction is set to zero and pressure becomes a function of x only. 

With the velocity component in the streamwise direction at the edge of 
the boundary layer, u g , specified as a boundary condition, dp/dx can be 
evaluated from the inviscid flow equation given by 

, du 

dp _ e 

dx ^e U e dx (2.27) 

The boundary- layer approximation provides very important 
mathematical advantages: first, the equations become parabolic instead 
of elliptic so that the streamwise direction becomes the marching 
direction and the numerical solution procedure becomes much simpler; 
second, pressure can be impressed upon the boundary layer as a known 
variable; third, boundary conditions are reduced considerably, notably 
for v and in the x coordinate direction. 

Note that in the case of turbulent flow, new unknown terms are 
introduced in the boundary- layer equations due to the Reynolds 
decomposition and averaging process. These terms - u'v' and - H'v', 
representing the apparent turbulent shear stress and heat flux 
quantities, must be modeled using empirical information to close the 
system of equations. The study of turbulence modeling will be discussed 


in Section II. E. 
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2. Boundary conditions 

Appropriate boundary conditions are required to solve these 
governing equations for viscous flow regions. In a standard procedure, 
initial profiles for u and H are required at the starting plane as well 
as values of u, v and H at the wall and values of u and H at the 
boundary layer edge. 

This standard procedure, which is referred as a direct method, 
becomes singular at the point of separation (Goldstein, 1948; Brown and 
Stewartson, 1969). When the pressure gradient is fixed near separation, 
the normal component of velocity, v, increases toward infinity at the 
point of separation. In a finite-difference solution, the magnitude of 
v at the point of separation is finite with a finite streamwise step 
size but will increase as the step size is reduced. As a result, the 
solution will not be unique and will usually fail to converge. This 
classical separation singularity, which is purely mathematical, can be 
avoided by the use of an inverse method which was suggested by Catherall 
and Mangier (1966). It has been used successfully by several 
investigators in numerical calculations (Klineberg and Steger, 1974; 
Williams, 1975; Carter, 1978; Kwon and Pletcher, 1979). In the inverse 
method, a regular solution is obtained through separated flow regions by 
prescribing the displacement thickness or skin friction in place of 
pressure gradient or u g and the pressure gradient is determined as a 
part of solution. 

The initial profiles for velocity and enthalpy at the stagnation 
point are provided automatically by similarity solutions to the 
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transformed equations that will be presented in a later section. The 
typical initial guesses required to obtain these solutions were obtained 
from freestream conditions, i.e., for y > -t, 

u(0,y) = u v(0,y) = 0 H(0,y) = H (2.28) 

At the solid boundary, the no slip condition was used for the 
velocity components: 

u(x,-t) = v(x,-t) = 0 (2.29) 

For the total enthalpy, either the wall value or wall heat flux was 
specified, i.e., 

H (x>-t) = H w (x) or = q w (x) (2.30) 

y t 

In the direct method, values of the velocity and total enthalpy 
were prescribed at the outer edge of the boundary layer as 

as y ■* ~, u(x,y) -*■ u (x) H(x,y) -*• H (x) (2.31) 

where the subscript e refers to conditions at the edge of the boundary 
layer. The pressure gradient was also determined from the specified 
u g (x) as given by Equation (2.27). 

The inverse procedure was used for the region interacted with the 
inviscid flow solution by specifying the displacement thickness which is 
defined as 



(2.32) 
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for axisymmetric compressible flows instead of the value of u at the 
edge of the boundary layer. The total enthalpy was specified as in the 
direct method. 

Although the inverse method was developed mainly for separated 
flows, the method can also be used for attached flows. However, in the 
present study, attached flows were generally computed by the direct 
method. From now on, for simplicity, the bars and tildes will be 
omitted from the single mean quantities. 

E. Turbulence Modeling 

In a previous section the need for turbulence modeling was pointed 
out. Unfortunately, a universal turbulence model has not been developed 
to date and it seems very unlikely that one will be developed soon. 
Therefore, it might be better to find turbulence models which have 
reasonable accuracy over limited ranges of flow conditions. With such a 
purpose, it is essential to understand the basic aspects of the 
structure of turbulent flow before proceeding with the implementation of 
turbulence models. 

1. The structure of the turbulent boundary layer 

As discussed in a previous section, the structure of turbulent flow 
appears to remain almost unchanged for Mach numbers up through 5. It 
is, therefore, generally sufficient to include compressibility effects 
implicitly through the mean density variation in turbulence models which 
work reasonably well for incompressible flows. 


c c 




The inner and outer regions of an incompressible, constant property 
turbulent boundary layer along a solid wall generally have quite 
different characteristics. However, they are strongly coupled by the 
shear stress profile and general diffusivity of the turbulence. The 
inner region comprises only a small fraction of a whole boundary layer 
in terms of thickness, but its influence over the entire boundary layer 
is significant. Klebanoff and Diehl (1952) experimentally observed that 
the inner region is generally insensitive to flow conditions far away 
from the wall and to the upstream conditions and that the mean velocity 
distribution is strongly dependent on the local conditions such as the 
wall shear, T^, density, p, viscosity, y, and the distance y from the 
wall. This suggests that the inner region is frequently in a state of 
near local equilibrium according to Bradshaw (1972). In such a case, 
the mean flow motion can be described by a rather simple expression 
known as the law of the wall, 

+ - , 

u = f x (y ) (2.33) 

by using the wall coordinates defined as 


+ u 
u = — 
u 


y + = yy 


(2.34) 


where u^ , which is called the friction velocity, is given by 


u 

t 



(2.35) 


The inner region can be further divided into three layers: a viscous 
sublayer, a buffer region, and a fully turbulent region based on the 
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relative magnitude of the viscous and turbulent shear stress. 

The outer region of the turbulent boundary layer occupies most 
(- 80%) of the boundary- layer thickness. The mean velocity distribution 
in this region is generally known to be described by the velocity-defect 
law: 


u - u 
e 


u 

T 


f 2 ( 6 > 


(2.36) 


where 6 is the boundary- layer thickness. Unlike the inner region, the 
function f is strongly affected by streamwise pressure gradient. 

Clauser (1956) discovered that a similarity velocity profile in the 
outer region can be obtained by choosing proper scaling variables. 

The preceding argument is generally applicable to fully attached 
turbulent boundary layers. When the turbulent flow separates under a 
strong adverse pressure gradient, the flow is often found to be 
unsteady, sometimes randomly and sometimes in a quasi-periodic sense. 
Simpson et al . (1977) observed that the qualitative turbulence structure 

upstream of the separation is not significantly different from that in 
flows with zero pressure gradient. As separation is approached, the 
flow is gradually influenced by the large scale outer flow and finally 
the motion near the wall is governed by the large eddy motion downstream 
of the separation point. These turbulent fluctuations in the separated 
region are of unusually large magnitude compared to the mean velocities. 
The law of the wall for the mean velocity profile and the local 
equilibrium argument appear not to be valid for separated flows 
according to Simpson (1979). 
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In the wall boundary layer developing downstream of the 
reattachment point, the mixing layer and the new wall shear layer 
intersect and result in a complicated turbulence structure. The 
incoming mixing layer in the outer part of the boundary layer is 
believed to carry characteristics of the separated region far downstream 
of reattachment. It has been also observed that it takes quite a long 
distance for flows to return to the structure of the ordinary turbulent 
boundary layer (Bradshaw and Wong, 1972). After the reattachment point, 
the turbulent energy is usually decreasing continuously but the exact 
reason for this phenomenon is not known. 

The above discussion was mainly focused on incompressible flows. 

As mentioned before, there is evidence that the basic structure of 
turbulence is not altered significantly by moderate density or 
temperature fluctuations, which suggests that interaction between the 
velocity and temperature fluctuations is probably not strong even in 
flow of moderate Mach numbers (Bradshaw, 1977). Thus for compressible 
turbulent shear flows, the main coupling between the governing equations 
occurs through the mean density variation only. 

As was the case for incompressible flows, the mean velocity 
distribution in the inner region can be observed to be somewhat similar 
by using wall coordinates. Since density and viscosity vary with the 
position, their values at the wall are used to define the friction 
velocity and wall coordinates as follows: 


+ 

u = 


u 

u 


, yu p 
v + = :.. t . w 

y 

w 


x i 

\ - <r ) 

r w 


(2.37) 
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It is, however, not easy to obtain a simple law of the wall expression 
for the mean velocity in the inner region in the compressible case. 
Analyses of experimental data suggest that law of the wall for 
compressible flows is significantly affected by compressibility, i.e., 
effect of high Mach number, pressure gradient and heat transfer at the 
wall. The usual law of the wall for compressible flows then becomes 

u+ = f 3 (y + » P + ) (2.38) 


where 




^c T p C u 
W W pw T 


u 

M = — 
t a 

w 


, v u du 
+ _ e e e 

13 u 3 dx 

T 


(2.39) 


and and a w are the outward heat transfer rate and speed of sound at 
the wall (Bradshaw, 1977). Generally, the effects of high Mach number, 
a hot wall (negative & c )» and a favorable pressure gradient (negative 
p + ) tend to drive the velocity data u + (y + ) down below the incompressible 
logarithmic law (White, 1974). There have been numerous suggestions for 
the law of the wall function, f (Van Driest, 1951; Maise and McDonald, 
1968; Bertram, 1968; Baldwin and MacCormack, 1976; Viegas and Rubesin, 
1983). These will not be given here and such a function was not used in 
the present calculations. 


Even though the basic structure of turbulence remains unchanged for 
flows with a moderate Mach number, the flow pattern might become quite 
complex due to the intersection of a shock wave with a boundary layer in 
transonic flows. In the inner part of the boundary layer near the wall 
where the flow is subsonic, disturbances created by the impingement of 


the shock wave cannot be discontinuous and thus are partly propagated 
upstream. This causes the streamlines upstream in the subsonic region 
to diverge and the increase of the thickness results in compression 
waves in the outer supersonic region. These intersect with the shock 
wave and cause it to bend forward. The pressure rise along the wall is 
still steep but continuous and takes places over a distance of the order 
of the boundary- layer thickness. If the pressure rise is sufficiently 
large , then the boundary layer separates with a more complex flow 
pattern forming a separation bubble and increasing the streamwise spread 
of the pressure rise. The complicated structure is primarily due to the 
mixing of supersonic and subsonic regions. 

During the interaction with the shock wave, relatively large shear 
stress gradients normal to the wall build up in the inner region and are 
of crucial importance, dragging low energy flow downstream into the 
region of high pressure. In the outer supersonic layer, the flow is 
more nearly inviscid and is usually described in terms of wave pattern. 
The most distinguishing features of this pattern are the refraction of 
all waves by the rotational supersonic layer and the reflection of all 
incoming waves when they reach the sonic line. The reflection at the 
sonic line is compatible with the behavior of the subsonic region 
(Green, 1970). 

If the overall pressure rise associated with the shock wave is not 
large enough to cause separation, the streamwise extent of the 
interaction region is then typically two or three times the thickness of 
the undisturbed boundary layer. In this case, since the upstream 
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propagation of pressure disturbances is small and the boundary layer 
remains thin, most of the outgoing compression waves coalesce with the 
outgoing random shock. Therefore, the effect of the viscous layer upon 
the pattern of the shock wave is relatively small. 

When separation occurs, the flow pattern becomes quite complicated. 
As shown in Figure 1, separation generally results in bifurcation of the 
shock wave into a leading and rear shock, generating a vortex sheet with 
a rapid entropy rise. A separation bubble of slow recirculating flow 
occurs at the foot of the shock wave in the lower portion of the viscous 
layer. A strong transonic interaction may involve the additional 
complication of a tongue of supersonic flow downstream of the shock. 

This supersonic region embedded in the subsonic region outside the 
boundary layer is believed to interfere with the normal process of 
reattachment. As a result, the length of the separated region becomes 
more sensitive to changes in the overall pressure rise. This appearance 
of an embedded supersonic zone is consequently associated with the 
beginning of the phase in which the interaction increasingly affects the 
pressure at the trailing edge. A more detailed description of the 
behavior of a transonic turbulent boundary layer subject to strong 
interaction with shock waves can be found in the articles by Green 
( 1968 , 1970 ). 

One of the major difficulties in predicting such a flow is 
uncertainty about the development of the turbulent shear stress in the 
region where pressure changes rapidly. This problem becomes serious 
when there is a severe adverse pressure gradient downstream of the 
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shock. In this case, the effect of uncertainties in predicting shear 
stress changes through the shock may be significantly magnified in the 
region of rising pressure downstream. Hence, as shock strength 
increases, it becomes more desirable to use a method which takes account 
of history, i.e., which includes an independent relation for the rate of 
change of shear stress. However, present understanding of a strongly 
interacted boundary layer is mainly qualitative and more detailed 
experimental studies are needed. Thus, present turbulence models are 
not complete and must be refined through a series of comparisons with 
experiments . 

2. Turbulence modeling 

The apparent turbulent stress and heat flux appearing in the 
Reynolds -averaged boundary- layer equations must be specified in order to 
predict the mean velocity and temperature distributions across the 
boundary layer. While relationships among these quantities can be 
developed into the form of transport equations from the basic 
conservation laws, more unknown quantities are introduced. These 
quantities must be evaluated based on empirical hypotheses. 

One of the simplest modeling strategies follows the proposal of 
Boussinesq (1877). Boussinesq assumed that turbulent shear stresses are 
related to the rate of mean strain via a turbulent viscosity defined by 


i i 3u 

■pu v = y — 
t 3y 


(2.40) 


By analogy with the kinetic theory of gases, which provides an accurate 
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theory for the molecular viscosity, the turbulent viscosity, , is 
assumed to be the product of a velocity and a length, 


y t = pv t* (2.41) 

where v^ and i are characteristic velocity and length scales of the 
turbulence, respectively. Therefore, the main task in representing the 
turbulent viscosity is to find appropriate expression for v^ and l. 

In a like manner, the apparent turbulent heat flux is related to 
the turbulent viscosity and the mean flow variables through Boussinesq*s 
turbulent conductivity concept. Using Equation (2.7), the apparent 
turbulent heat flux in the boundary layer is approximated as 

pH T v ! = ph T v T + puu f v f (2.42) 


Using Boussinesq's turbulent conductivity concept and the turbulent 
Prandtl number, 


-phV 


y t 3h 
Pr t 3y 


(2.43) 


Substituting Equations (2.40) and (2.43) into Equation (2.42), the 
apparent turbulent heat flux is given as 


77T-, V t 8H ^ ,, 1 . 3u 

-pH v = — — — + (1 - - — )y u— 
r Pr 9y Pr t t 


(2.44) 


Experimental data suggest that Pr is a well behaved function across the 
flow, and reasonably good predictions have been achieved even with a 
constant value of for example, 0.9. Note that yi^ and Pr^ are not 

the physical properties but vary with local flow conditions and 


geometry. Although the Boussinesq assumption is not in complete 
agreement with all available experimental evidence, it appears to be a 
reasonably good approximation in many engineering flow circumstances. 

The models that utilize the Boussinesq assumption are further 
classified by the number of supplementary partial differential equations 
which are used to obtain the modeling parameters like v fc , Z or 
itself. An ordinary differential equation is usually counted as a half 
equation and algebraic equations are counted as zero. 

The algebraic or zero equation model, which is based on Prandtl’s 
mixing length concept (1926), is the simplest and most popular 
turbulence model among those utilizing the Boussinesq assumption. In 
this model, the characteristic scales and the turbulent Prandtl number 
are given by simple algebraic equations related to the motion of the 
mean flow. Despite its simplicity, it has proven effective in 
predicting relatively simple flows. In order to make it more accurate 
for complicated flow cases, there have been numerous attempts to modify 
the algebraic equations with semi-empirical relations (Van Driest, 1951; 
Patankar and Spalding, 1967; Cebeci and Smith, 1974; Deiwert, 1976; 
Baldwin and Lomax, 1978; see, e.g. , Anderson et al . , 1984). The major 
objection to algebraic models is that they are based on a local 
equilibrium assumption, i.e., the turbulent viscosity is evaluated only 
in terms of local flow parameters and upstream effects are not 
considered . 

A one-half equation model is considered as the least complex method 
which can approximately account for the nonequilibrium effects. One of 
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the modeling parameters or turbulent viscosity itself is controlled by 
the solution of an ordinary differential equation which can usually be 
derived from the more general transport equation for the parameters by 
neglecting variations in one coordinate direction while the other 
parameters, if any, are governed by algebraic relations. Most of these 
models employ an ordinary differential equation for the length scale 
(McDonald and Kreskovsky, 1974; Chan, 1972; Pletcher, 1978) or turbulent 
viscosity (Shang and Hankey, 1975; Reyhner, 1968). 

Obviously, the next step is to use a full partial differential 
equation for the modeling parameters. The modeling which has become 
popular in not only academic research but also engineering applications 
is a two-equation model. One of the most frequently used two-equation 
model is the k -z model first proposed by Harlow and Nakayama (1967) and 
developed further by many others (Jones and Launder, 1972; Ng and 
Spalding, 1972). Also numerous other two-equation models have been 
suggested including a k-u> 2 model developed by Wilcox and Rubesin (1980) 
especially for compressible flows. Most of these two-equation models 
employ a modeled form of the turbulent kinetic energy equation but use a 
different dependent variable for the second model transport equation 
from which the length scale is determined. The disadvantage of two- 
equation models is the need to make assumptions in evaluating the third- 
order turbulent correlations in the transport equations. 

A Reynolds stress equation model which does not utilize the 
Boussinesq assumption about turbulent viscosity was pioneered by Rotta 
(1951) and has been enhanced by Hanjalic and Launder (1972) and many 


others. Recently, a large eddy simulation approach which is not based 
on the Reynolds equations was developed by Deardorff (1970). However, 
these last two types of models must be refined and tested further before 
they can be used for engineering predictions. 

Transonic turbulent flow is believed to be affected strongly by 
upstream history effects in the neighborhood of a shock wave. For such 
flows, two-equation models have been used extensively because they were 
believed to have a better chance of predicting overall flow structure 
correctly than local equilibrium models. Surprisingly, most of 
calculated results of various two-equation models (Coakley and Viegas , 
1977; Viegas and Coakley, 1978; Viegas and Horstman, 1979) are not 
necessarily superior to those of lower-order models. Meanwhile, several 
one-half equation models have been developed mainly in order to account 
for nonequilibrium effects and their predictions generally provide a 
definite improvement over those from algebraic models. In the present 
study, the zero equation and the one-half equation models were used and 
they are described in detail in the following sections. 


3. Algebraic model 

Prandtl (1926) proposed the following mixing- length formulation: 




t 



3u 

9y 


(2.45) 


where £ and £|3u/3y| can be thought of as the characteristic length and 
velocity scales of turbulence, respectively. Evaluation of £ in this 
model varies with the type of flow. For flows along a solid surface, 
the boundary layer is usually divided into inner and outer regions for 
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the evaluation of £. 

In the inner region close to a solid wall, the mixing length is 
believed to be proportional to the transverse distance from the wall as 

£ . = icy (2.46) 

1 

where k is an empirical constant known as the von Karman constant. A 
value of k of about 0.41 provides good agreement with experimental data 
taken in simple flows. This distribution of £^ matches correctly with 
the fully turbulent layer of the inner region excluding the viscous 
sublayer and buffer layer close to a wall. The expression for mixing 
length become applicable over the whole inner region by using an 
empirically determined damping function D, proposed by Van Driest 
(1956). 


£^ = icDy 


(2.47) 


where 

D = 1 


4 - 

exp(- ^-) 
A 


(2.48) 


and y + is given by Equation (2.39). A value of 26 is generally used for 
the damping constant A + . The expression given by Equation (2.48) is 
valid only for flows with a negligible pressure gradient. When the 
pressure gradient is sufficiently adverse as to cause separation, 
becomes zero, then D given by Equation (2.48) will become zero so that 
will be zero in the inner region. Such a problem can be overcome by 
defining D based on the maximum velocity gradient rather than the wall 
value as suggested by Pletcher (1978) and Carter and Wornom (1975). 
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Then the damping function can be written as 


D = 1 - exp [ - ( 


-fB. 


3u 

3y 


* y 

> 


(2.49) 


' max A 

There have been numerous variations suggested for the damping 
function in order to account for effects of property variations, 
pressure gradient, blowing, transverse curvature and separation (Cebeci, 
1970, 1971, 1973; Kays et al., 1970; Pletcher, 1976). Following Cebeci 
and Smith (1974) the damping function used in the present study is 


D = 1 - expl-Cj 


9u 

3y 


>* *] 

max A 


(2.50) 


where 


N = [1 - n. 8 A)(fi-) 2 pV 
p w 


(2.51) 


In the outer region of a wall boundary layer, the mixing length is 
often taken as proportional to the boundary- layer thickness. Another 
common method is to use the kinematic displacement thickness instead of 
boundary- layer thickness and the velocity at the outer edge of the 
boundary layer as the length and velocity scales, respectively: 


V. = apu 5. 
to r e k 


(2.52) 


where a is the Clauser constant and is observed to be about 0.0168 and 
6, is the kinematic displacement thickness defined by 


\ = ' r (1 - ib 

0 o e 


(2.53) 


6o 


Cebeci and Smith (1974) also recommend a modification to include 
the effect of low Reynolds number, based on the observation by Coles 
(1962) as 

1.55 


a = 0.0168 


1 + ir 


(2.54) 


where 


tt = 0.55 [1 - exp(- 0.243z 2 - 0.298z)] 


(2.55) 


and 


Re 


z = 


Ok 


425 


- 1 


(2.56) 


Here Reg^ is the Reynolds number based on the kinematic momentum 
thickness defined as 


Re ek = 


p U 0 

r w e k 


w 


where 0, is the kinematic momentum thickness defined by 
k 


0 t r u u . , 

k r u u ' J 

0 o e e 


(2.57) 


(2.58) 


The effect of this modification becomes negligible if Re^ is greater 
than 5000. 

As the freestream is approached near the edge of the boundary 
layer, the turbulence becomes intermittent. Following Klebanoff ! s 
(1954) observations, there are intervals of time when the flow is not 
turbulent near the outer edge of the boundary layer and these intervals 
become longer as the distance from the wall increases. Including this 
intermittency factor to the expression for turbulent viscosity for the 
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outer region. Equation (2.56) becomes 


u = apu 6, X . 
to r e k i 


where intermittency factor 2f\ is given by 


6 -1 


= [1 + 5.5(|) ] 


(2.59) 


(2.60) 


and a is given by Equation (2.54). Note that the effect of the 
transverse curvature was neglected since the boundary- layer thickness is 
assumed small compared with the radius of the body. Deiwert (1976) used 
the furthest point across the boundary layer where the velocity is zero 

in place of usual y = 0 for the lower basis needed in evaluating 6^ and 
in order to avoid an unrealistically large length scale in the 
reversed flow region. The effect of this modification proved to be 
insignificant and was not used in the present study. 

The above algebraic two-layer model will be referred to as the 
Cebeci-Smith model hereafter. This Cebeci-Smith model is known to 
predict rather poorly in separated regions and was not proposed for use 
in such flows by its originators. This poor prediction is mainly due to 
the assumption of local equilibrium which implies the neglect of the 
effects of diffusion and convection of turbulence scales and assumes a 
balance between the production and dissipation of the turbulent kinetic 
energy. The sole reason for its use was for reference purposes, simply 
to indicate the level of performance to be expected from algebraic 
models which have been perfected for equilibrium or near -equilibrium 


turbulent flows. 


70 


4. Nonequilibrium turbulence model 

The aforementioned algebraic model is based on an assumed 
equilibrium between the mean motion of the flow and its local 
turbulence. When the flow changes rapidly in the streamwise direction, 
the turbulence may lag the mean motion. To account for this, models 
have been developed to utilize transport equations for the turbulence 
itself. The simpler of such methods use transport equations for 
variables used in defining a turbulent viscosity, which then can be out 
of equilibrium with the mean motion. It should be noted, however, that 
the Reynolds stresses are tied to the mean motion through the Boussinesq 
approximation and react immediately to changes in the mean motion, even 
though the effect of this reaction is influenced by the extent of the 
lag of the turbulent viscosity. 

Such models typically employ partial differential equations to take 
into account the effect of diffusion or convection of turbulence length 
or velocity scales. Typical of such models is the two-equation model 
(e.g., k-£ model of Jones and Launder (1972) and k-w 2 model of Wilcox 
and Rubesin (1980)). Considering the added numerical complexity and 
computational effort required to solve additional partial differential 
equations, the improvement in predictive accuracy over algebraic models 
has not been encouraging. On the other hand, the effect of flow history 
on the turbulent viscosity can be approximated with an ordinary 
differential equation, a so called one-half equation model, and such 
models have proven effective within a limited range of applicability. 
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One-half equation models are generally divided into two different 
types. The first type is purely an empirical relaxation or lag model. 
The models of this kind are actually equivalent to one-dimensional 
versions of transport equations for the quantities concerned except that 
these transport equations are not generally derivable from the Navier- 
Stokes equations. Included in this type are a relaxation form of an ODE 
for itself proposed by Shang and Hankey (1975) and an ODE for length 
scale in the outer region suggested by Pletcher (1978). The other type 
utilizes a reduced form of the transport equation for turbulent kinetic 
energy to define one of the turbulent characteristic scales (McDonald 
and Kreskovsky, 1974; Chan, 1972). 

Very recently, Johnson and King (1985) suggested a new one-half 
equation model in which the velocity scale is governed by an ODE 
derivable from the turbulent kinetic energy equation. This model will 
be referred to as the Johnson-King model and is described below. 

The turbulent viscosity is assumed to have a functional form of 


V t = y to [1 ‘ exp(_ 7^1 (2 - 61) 

H to 

where y^ and y can be thought of as describing the turbulent 
viscosity in the inner and outer parts of the boundary layer. Equation 
(2.61) provides a smooth implicit blending of turbulent viscosities in 
two regions instead of an explicit change at the boundary of two 


regions . 

The inner turbulent viscosity is given by 
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2 4 

y ti = pU icy(-u'v' m ) 2 (2.62) 

where the damping function D is the same as given by Equation (2.50) 
except that A + is assigned a value of 15 instead of 26. This adjustment 

4 

was suggested by Johnson and King (1985) to provide y -dependency 
3 

instead of y -dependency of in the near wall region. Compared with 
the algebraic Cebeci-Smith model, the major difference is that the 
velocity scale is based on (-u'v' )^, which is provided by the solution 
of an ODE. Perry and Schofield (1973) suggested that (-u'v'^) 2 provided 
a much better velocity-defect correlation than did the friction velocity 
u^. for flows with adverse pressure gradients. 

The outer value of the turbulent viscosity y is given by 

\o = ° (x)au e 6 kh (2 - 63) 

where o(x) is an unknown parameter which varies with streamwise 
distance. This form of the outer turbulent viscosity is identical with 
the Cebeci-Smith model except for the appearance of o(x) which becomes 
necessary to make the turbulent viscosity consistent with the local 
flow. In other words, values of o(x) are determined so that -u'v* 

m 

resulting from use of y^ and the mean velocity profile coincides with 

the value given by the ODE. In this way, the nonequilibrium effect 

expressed through -u'v' can be felt in the outer region as well as in 

m 

the inner region. The details are discussed later. 

The streamwise distribution of -u f v* is determined from the 

m 

solution of an ordinary differential equation which is derived from the 
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turbulent kinetic energy equation using assumptions similar to those 
used by Bradshaw et al. (1967). The equation for the transport of 
turbulent kinetic energy for a steady compressible boundary layer is 
given by 


9k 9k r i 9 r 9 n , *2\i 

pU 3^ + PV ^ = _pU V ^ + 3^ [y 3^ (k + V )] 


— (pkv' + pv') - pe d 


(2.64) 


where k is the turbulent kinetic energy and is defined as 


k = \ (u’.up 


(2.65) 


and s , is the rate of dissipation of the turbulent kinetic energy and 
d 

can be approximated as 


9u! 2 

Hr— i) 

P^x/ 

J 


( 2 . 66 ) 


Assuming the path of the maximum turbulent kinetic energy is 
continuous outside the viscous sublayer and is nearly coincident with 
the x coordinate, Equation (2.64) can be specialized for the maximum 
turbulent kinetic energy k^ and reduced approximately to the following 
form: 


dK m -r~. 3u| 1 a , —< .L — ; n 

u j — = -u v — - - — (pkv + pv ) - £ , 

m dx m 9y| m p 3y r r d,m 


(2.67) 


Here, subscript m denotes that the quantity is evaluated where k is 
maximum in the y direction across the boundary layer. In the above 


equation, the additional assumption has been made that -u v becomes 


maximum where k is maximum. 
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Following the analysis of Bradshaw et al. (1967), additional 
parameters are defined as 


-u' v' 


3/2 


a l = 


m 


m 


L = 
m 


(-uV ) 

m 

8 d,m 


( 2 . 68 ) 


Here L corresponds to a dissipation length scale. To simplify the 
m 

above equation, a^ is assumed to be constant. Experimental data suggest 
that a^ varies between 0.2 and 0.3 and it is reasonable to assume a^ to 
be a constant. A value of 0.25 for a^ was suggested by Johnson and King 
(1985). Different values of a^ were also used in the present 
calculations, and the effect of varying a^ will be discussed in a later 
section. Using Equation (2.68), Equation (2.67) can be rearranged to 


give 


a (-u v ) 
1 m 


^(-u’vV = L u 

m m 


[L S - C--V. >* - ^ » 1 


3u 
"m 3y 


m 


u m J 
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(2.69) 


where D represents the diffusion term in Equation (2.67). 
m 

The first term on the right hand side of Equation (2.69), 

L (3u/3y)| , can be interpreted as the square root of the turbulent 

m m 

shear stress that might result when convection and diffusion effects are 
negligibly small. Thus, this term is replaced as 


3u I 


= (-U v ) 
m 3y m,eq 

m n 


(2.70) 


This quantity is determined by the equilibrium turbulent viscosity 


distribution given by 
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V = V. [1 - exp( 3)1 

t,eq to,eq l y Ji 

to.eq 


(2.71-a) 


y ti,eq 

y to,eq 


= pD Ky(-u'v' ) : 
m,eq 


= au 6. % . 
e k i 


(2.71-b) 

(2.71-c) 


The dissipation length scale L is determined as 

m 

L m = 0,4 y m for y m ~ °’ 225 6 (2.72-a) 

L m = 0-09 6 for y m > °' 225 6 (2.72-b) 

The turbulent diffusion term must also be modeled in order to 
solve Equation (2.69). This turbulent diffusion term seems to be 
important in the region where the flow is recovering toward an 
equilibrium state. Johnson and King (1985) proposed the following form 
of D m based on the bulk convection hypothesis of Townsend (1976): 

3/2 

C , (-u V ) . 

n _ dif nr_ ,, , 

D m a.6[0.7-(y/6) ] I 1 ‘ o(x) I (2.73) 

where is a modeling constant and its suggested value is 0.5. 

In order to solve Equation (2.69), it is convenient to define new 
variables : 


8 = (-uV ) 

m 


2 0 =(-u'v' ) 2 

eq m,eq 


(2.74) 
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Finally, Equation (2.69) can be rewritten in terms of new variables as 
follows : 


d9 _ a l f 8 C dif L m 

dx 2u L U ' 9 a.6[0.7-(y/6) ] 

mm eq 1 m 


1 - o (x) 2 | } 


(2.75) 


Equation (2.75) is nonlinear because of the appearance of local 
parameters as coefficients of 0. Equation (2.69) is linearized by 
simply using values of those quantities from the previous x station 
(lagging). Also, in Equations (2.62) and (2.71-b), the values of 


* * )2 


(-u’v’ ) 2 and (-u r v T ) 2 required to determine u . and u . are 
v m m, eq n ti ti,eq 

obtained from the previous x station. 

Initial conditions for 0, 0 and o(x) must be provided before 

Equation (2.75) can be solved. The flow upstream at the point where 

this nonequilibrium model starts can be assumed to be in an equilibrium 

state so that u = y . Therefore, the initial condition is given as 
t t , eq 


at x = x , 
o’ 


u ! v* = u ! v ! (i.e. , 0 = 0 ) 

m m , eq eq 


(2.76) 


o(x) = 1.0 


where -u’v' can be easily provided by any equilibrium turbulent 

m , eq 

viscosity model such as the Cebeci-Smith model. The starting point, x , 
for this closure should not, obviously, be where the Reynolds stress is 
expected to change rapidly. Therefore, all calculations were initiated 
using the Cebeci-Smith model from the leading edge and a switch was made 
to the Johnson-King model at the beginning of the viscous - inviscid 


interaction region . 


/ / 

The variable o(x) provides a close tie between the turbulent 
viscosity across the whole boundary layer and the ordinary differential 
equation for the maximum Reynolds shear stress. The appropriate value 
of o (x) can be obtained through an iterative process by requiring that 
-uV m obtained from the velocity profile and turbulent viscosity agrees 
with -uV m determined from Equation (2.75). To evaluate y , 
specifically an initial guess for o(x) is necessary to begin the 

iteration process. Normally, the value of o(x) from the previous 
station was used. Then the intermediate distribution of y noted by 

y t , is determined by Equation (2.61) and the distribution of Reynolds 

“V 

stresses, -u v , based on y^_ is calculated as 



(2.77) 


If -uV m obtained by Equation (2.77) is different from -uV 

m 

determined from Equation (2.75), the value of o(x) must be modified. 

The velocity gradient at the point where -VV* becomes maximum is then 
given by 


3u 

8y 


m 


t T 

~pu v 

m 

t y m 


(2.78) 


The turbulent viscosity which gives the same maximum Reynolds stress 
obtained by Equation (2.75) should be 



(2.79) 
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The new value of o(x) can then be determined by letting 


U* Jo) = VL m ( 2 - 80 ) 

t y m t y m 

Since the expression for y is nonlinear in terms of o(x), the solution 
to Equation (2.80) has to be determined by a numerical method. In this 
study, the simple Newton-Raphson method was applied as 


k+1 k 

Ti — y 

to,m to,m 


U. m “ V,. m 
t , m t , m 


where superscript k is the iteration level for the Newton-Raphson 

k+1 

method. The converged value of u becomes (u ) and new 

to,m to,m new 

value of o(x) is given by 

(y ) 

to,m new /0 c 

o = 0 - , t r (2.£ 

new old (y m ) n , 
to,m old 

Since y in Equation (2.79) is still dependent on y , i.e., the o] 
t ,m t ,m 

value of o(x), this procedure of determining o(x) needs to be 
iteratively repeated . 

The solution procedure for this closure can be summarized as 
follows : 

1) For x < x q , the flow is assumed to be in equilibrium so that the 
turbulent viscosity is provided by the Cebeci-Smith model. 


2) At x = x^, initial conditions are given by Equation (2.76) and 

-u’v' is determined by the Cebeci-Smith model, 

m , eq 

3) For x > x , y and -u ! v* are determined by Equation (2.71). 

o t,eq m , eq 

Also, o(x) is lagged and -uV is calculated by Equation (2.75). 

m 


7 a 

/ v 

* -v 

4) y is calculated by Equation (2.61) and -u'v' is determined bv 

^ m J 

JU 

Equation (2.77) using y t and the velocity profile. 

5) If -u'v' A -u'v' , the new a(x) is calculated by Equation (2.82) 
and step (4) is repeated. Otherwise the procedure moves to step (6). 

“V 

6) The turbulent viscosity is determined by y^_ = y^_ 

This Johnson-King model has never been tested for a wide range of 
flows, but appears to be very effective in predicting both pressure 
driven and shock induced separated flows in the present study. However, 
this model still needs improvement in order to provide a better 
prediction and this will be discussed in a later section. 

Before closing the discussion of turbulence modeling, it should be 
noted that for turbulent flows, transition from laminar to turbulent 
flow was assumed to occur at a fixed point near the leading edge and no 
particular formula for the transition point was used. Also, the 
apparent turbulent heat flux is modeled by Equation (2.44) with a 
constant value of the turbulent Prandtl number equal to 0.9. 

F. Mathematical Model 

Although the boundary- layer equations given by Equations 
(2 . 24) - (2 . 26) are believed to describe the compressible turbulent 
boundary layer fairly well, they are not in the best form to be solved 
numerically. The additional modifications and transformations used 
prior to obtaining the numerical solution are described next. 


c ~ ^ 
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1. Prandtl’s transposition theorem 

The shear- layer coordinate system has been adopted because of 
several advantages pointed out in Section II. C. The difficulty 
associated with the use of the shear- layer coordinate is that the body 
surface does not lie along a coordinate line. According to Werle and 
Verdon (1980), this difficulty can be avoided by using the Prandtl’s 
transposition theorem (1938). A detailed discussion of this theorem can 
be found in a article by Glauert (1957). In this transposition theorem, 
new variables defined by 

x - x y = y + t v = v + u~ (2.83) 

are introduced into the boundary- layer equations, where t is defined as 
the distance between the shear- layer coordinate and the body coordinate 
as shown in Figure 2. 

After application of this transposition theorem, the form of the 
boundary- layer equations does not change and the same boundary 
conditions are maintained except for the appearance of a caret (A) over 
x, y and v. Provided that the magnitude of t is only a fraction of the 
boundary- layer thickness, the boundary- layer approximation remains 
valid. The only differences are that the wall boundary conditions, 
given by Equations (2.29) and (2.30), are prescribed at y = 0 in place 
of y = -t and the velocity v is no longer the component of velocity in 
the y direction but becomes the component of velocity normal to lines of 
constant y. Therefore, the caret (a) will be dropped for x, y and v in 


the equations to follow. 
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Nondimens ionalization 

Variables in the boundary- layer equations are nondimensionalized as 


X = 


Y = ? 


Re 


R = r 


U = — 

U 


V = ^Re } 

u <*> 


I = 


u 


(2.84) 


P = — ^ 
P u 




~ y 
y = — 
y 




where Re^ is the Reynolds number based on the freestream value and L is 
the reference length of unity. Also note that the tildes (~) over the 
density and viscosities do not mean a mass-weighted average. 

Using these nondimens ional variables and Boussinesq's turbulent 
viscosity and turbulent conductivity concepts (see Equations (2.40) and 
(2.44)), the compressible boundary- layer equations for steady, 
axisymmetric flows, given by Equations (2 . 24) - (2 . 26) , become 
continuity 


fjj(pUR) + fy(pVR) = 0 


(2.85) 


momentum 


~ T t3U ~ tt 3U _ ~ n , 1 8 r ~ ~ N 9U, 

pu 3X " pv 3Y p e U e dX ^ R 3Y tR(v + y t ; 3Y ] 


(2.86) 


gg grg y 


~..3I ~.3I 

pU 5 x ' pV 3 Y 


I i_ (Rf L + 

r 3Y w Pr Pt 'ay 


+ R[(i - (i - p^)v t ]u|H) 


(2.87) 
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The corresponding boundary conditions, given by Equations 


(2.31), can 

be rewritten 

as 



at X = 0, 

U = 1.0 

V = 0 

i = i 

oo 

(2.88) 

at Y = 0, 

U = V = 0 



(2.89) 


I = I 

w 

or 

^ specified 

w 


as Y «>, 

u -*> U 

e 

I + I 

e 


(2.90) 


3. Coordinate transformation 

The boundary- layer equations in the form of Equations (2 . 85 ) - (2 . 87 ) 
are formidable nonlinear coupled partial differential equations even for 
a very simple geometry. Therefore, there have been a number of efforts 
to simplify the equations, especially in the pre-computer era. Among 
these earlier works, a number of transformations relating compressible 
boundary layers to equivalent incompressible flows can be found useful 
even today. 

The density variation can be accounted for by using a stream 
function, but the variable viscosity requires the introduction of a 
second similarity type variable. It was discovered by Illingworth 
(1950) that it was convenient to account for viscosity and density 
effects separately in the similarity variables. Lees (1956) combined 
the Illingworth transformation and Mangier transformation (1948) for 
axisymmetric flow and this transformation is often called the 
I llingworth-Levy-Lees-Dorodnitsyn transformation . This transformation 


was modified by including the transverse curvature terms following the 
lead of Probstein and Elliott (1956). Such a transformation was used by 
Cebeci and Smith (1974), Carter (1978, 1981) and Lee and Pletcher (1985) 
for the analysis of compressible turbulent flows. 

The stream function for compressible axisymmetric flow is defined 
as 


— - - 7?vr ^ - ~tu? 

ax “ pVR ay ~ pUR (2.91) 

Use of equation (2.91) replaces the continuity equation. The similarity 
type independent variables are defined as 


X 


5 = / py u R dx 

e e e o 


U Y 

T| = -tIt S pRdY 


8 ( 5 ) 


(2.92) 


0 0 

where the normalizing function, g(£), can be chosen arbitrarily. This 
transformation removes most of the effects of compressibility from the 
governing equations. Following Illingworth's work, when the direct 
method was used, g(£) is given by 


8 ( 5 ) = ( 2 ?) 


* 


(2.93) 


so that the singularity at £ = 0 can be removed. As a result, the 
calculation can be started without difficulty at a leading edge or at a 
stagnation point. When the inverse method was used, g(£) was taken to 
be equal to the displacement thickness (Carter, 1978). This choice 
tended to keep the number of grid points within the boundary layer 
nearly constant. 
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Defining new dependent variables. 



Equations (2 . 85 ) - (2 . 87) become 
continuity 


gF 


dtp 

3ti 


(2.94) 


(2.95) 


momentum 


2 3F dtp 3F _ 2 r 2 9_ ( V3F, 

* ^ " 8 &(G F } aT! [C e cl + V ^ 3 ti 


(2.96) 


energy 


2 3G dtp dG _ a_ r ^£ n U t Pr 3G 

8 *35 ’ S 3£ 3 ti 3n l Pr U v Pr t J 3Ti J 


+ C E - & + (1 - pj;»iri F S> 


(2.97) 


where the nondimens ional pressure gradient parameter, f), and the 
Chapman-Rubes in (1949) type parameter, C^, are defined as 



dM 
e 


c = -fiW-fB-) 

°fc p y l -R } 
e e o 


(2.98) 


and assuming is 


constant , 


the Eckert type parameter C £ is given by 


U 2 (Y-l) M 2 

e e 

C E = I = 1 + 0.5(y-l)M 2 
e e 


(2.99) 


Care must be taken in evaluating the convective terms in Equations 
(2.96) and (2.97) at the leading edge when $ = 0. The limiting value of 
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2 8 

§ terms is zero as £ - 0; however, the g^| ~ terms in 
Equations (2.96) and (2.97) has a finite limit as $ -*• 0. To obtain the 
most accurate numerical solution at the leading edge, a new dependent 
variable, W, defined as 


W = - 




( 2 . 100 ) 


was introduced into Equations (2.96) and (2.97) and the continuity 
equation, Equation (2.95), was reformulated from Equation (2.85) using W 
instead of rp. This procedure was used only at the leading edge. The 
stream function formulation given by Equations (2 . 95) - (2 . 97) was used 
thereafter. The details of solution procedure at the leading edge are 
quite similar to those used for the stream function formulation and will 
not be discussed further. 


This transformation was not used to solve the transport equation 
for the maximum Reynolds stress in the Johnson and King turbulence model 
because that transport equation is an ordinary differential equation and 
the transformation offered no advantage. 

The corresponding boundary conditions in the transformed 
coordinates for Equations (2 . 95) - (2 . 97) are given as 

At ti = 0, F = tf, = 0 (2.101) 


G = G 


w 


or 


9G 

3t)i 


specified 


w 


As T| 


F = G = 1.0 


( 2 . 102 ) 


The definition of the displacement thickness, Equation (2.32), gives the 
value of the stream function at the edge of the boundary layer as 
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T1 p 

l> e = g(5) / 6 (~)dn - m 


(2.103) 


where m is the nondimensional mass flux defined as 


tn = p U R T (2.104) 

r e e o 

and 6 is the nondimensional displacement thickness. 

Equations (2.96) and (2.97) are parabolic and can be solved in a 
forward marching manner in the streamwise, £, direction as long as 
F > 0. In regions of reversed flow, i.e., F < 0, the correct marching 
direction is the negative direction which is the true "streamwise" 
direction. With flow reversal, the solution to the complete boundary- 
layer equations can be advanced in the positive coordinate direction 
only iteratively and velocities in regions of reversed flow must be 
stored. In each iterative sweep, the correct difference representation 
for convective terms is provided through windward differencing which 
honors the appropriate marching direction. This method has been used 
successfully by many investigators (Klineberg and Steger, 1974; Carter 
and Wornom, 1975; Cebeci, 1976). 

On the other hand, a much simpler alternative to the multiple sweep 
procedure was proposed by Reyhner and Fliigge-Lotz (1968). According to 
their suggestion, which is known as the FLARE approximation, the 
streamwise convective term F80/8£ is replaced by C|F|90/8£ where C is 
zero or small positive constant if F < 0 and one if F > 0. This term 
has been observed to be very small in regions of reversed flow. The 
validity of the FLARE approximation in the case of small separation 


regions has been demonstrated in many calculations (Carter, 1974; 
Williams, 1975; Cebeci et al. , 1979; Kwon and Pletcher, 1981) and also 
confirmed by the experimental observations of Simpson et al. (1974). 
Recently, Davis and Carter (1984) showed that the FLARE approximation 
still gives results comparable to those obtained from the windward 
differencing scheme with reversed flow velocities even up to 28% of the 
edge velocity. Throughout the present study, the FLARE approximation 
was used for most of the calculations. The windward differencing scheme 
was also employed, but only for limited cases for purposes of 
comparison. 


G. Numerical Method 

This section describes the method used to solve the boundary™ layer 
equations developed in the previous section. The continuity and 
momentum equations were solved in a coupled manner using a fully 
implicit finite-difference scheme. This results in a block tridiagonal 
system with the blocks being square matrices of order 2. Kwon and 
Pletcher (1981) reported that coupling of the continuity and momentum 
equations eliminated the wiggles in the regions of separation which 
appeared in the solution when an uncoupled scheme was used. The energy 
equation and the additional transport equation used in the Johnson and 
King turbulence model were solved in an uncoupled manner. 

Following a discussion of the computational grid, the finite- 
difference discretization and solution procedures will be presented. 
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1. Computational grid 

A representative example of the finite-difference grid used for 
this study is shown in Figure 4. However, the mesh shown in Figure 4 is 
much coarser than those used in the actual calculation. The i index is 
used to specify the X-position. The point i = 1 corresponds to the 
leading edge. The j index is used to specify the Y-position with j = 1 

corresponding to the points at the wall and j = NJ corresponding to the 

points at the outer edge of the computational domain. The formulation 
of the grid can have a significant effect on the convergence and 
accuracy of the solution. The optimum method of grid formulation is 
usually achieved through trial and error and is different depending on 
the flow case. For laminar calculations, an equal grid spacing was used 
in the normal (ti) direction. For turbulent flows, wall functions were 
not used and the turbulence model was applied directly to the wall as 
pointed out earlier. Therefore, a variable r\ grid spacing was used so 

as to have more grid points concentrated near the wall and to avoid an 

excessive number of grid points. 

In this study, the grid scheme proposed by Cebeci and Smith (1974) 
was employed. This scheme maintains a constant ratio between two 
adjacent grid increments as 



where K is constant. In some cases, the accuracy of the solution was 
somewhat sensitive to the value of K used. A value between 1.05 and 
1.10 generally gave the most reliable results. In all turbulent cases, 
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the grid increment nearest the wall, Arj ^ (= r\^ - ti ) , was chosen so that 
the corresponding Ay + did not exceed 1.0, and a number of mesh points, 
usually 10, were placed within the viscous sublayer. 

Variable grid spacing was used in the streamwise (£) direction for 
both laminar and turbulent flows with most of the points concentrated in 
the region where the strongest pressure gradient was expected. A good 
example of this would be the concentration of points near an 
experimentally observed shock position. The simultaneous viscous- 
inviscid interaction method requires the same grid spacing in £ 
direction for both the boundary- layer equations and the potential 
equation. On the other hand, the two grids can be independent when the 
semi-inverse viscous - inviscid interaction method is used although the 
same grid was generally used in the present study. The grid spacing in 
the £ direction was generated using a stretching transformation similar 
to the one proposed by Roberts (1971). This will be discussed in detail 
in the inviscid analysis section. Outside of the interaction region, 
the streamwise grid spacing was chosen so as to be approximately 
proportional to the thickness of the boundary layer. A typical mesh 
size used 100 increments in the streamwise direction and 50 to 70 in the 
normal direction. 

2. Finite-difference representation 

The set of nondimens ional ized boundary- layer equations in the 
transformed coordinates discussed in Section II. F was solved over the 
region of interest using a fully implicit finite-difference method. In 
this section, the finite-difference scheme employed to solve the 
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boundary- layer equations and the equation for turbulence modeling is 
described. 

There have been many finite-difference methods proposed to solve 
the boundary- layer equations. These can be generally divided into 
explicit and implicit methods. The explicit method is simple but is 
often constrained by a stability condition, i.e., suffers from a severe 
limitation on marching step size. On the other hand, most implicit 
schemes in common usage are inherently stable. Among them, the Crank- 
Nicholson implicit scheme is second-order-accurate in both the 
streamwise and normal directions but barely satisfies the stability 
condition. This procedure has occasionally been found to be unstable 
for turbulent flows (Anderson et al., 1984). The fully implicit scheme 
used in the present study is first-order-accurate in the streamwise 
direction and second-order-accurate in the normal direction and is also 
unconditionally stable. When the fully implicit scheme is used, 
linearization and possible adjustments in the differencing in order to 
maintain diagonal dominance must be dealt with as in any other implicit 
scheme. The linearizing procedure is given in the next section. 

a. Continuity and momentum equations The finite-difference 
representation of the continuity and the momentum equations is given as 
follows . 
continuity 


g.„i . „i 

2 U j + *j-l J 




) 


(2 . 106) 
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M 1 = 
n j 




where 0 is a general dependent variable. A similar representation 
applies to the operators in the £ direction. 

The double arrow notation in Equation (2.108) on the £ difference 
operators means that the difference is always upwind, which is a 
backward difference when F > 0 and a forward difference when F < 0. The 
nondimens ional diffusion coefficient N. tl was evaluated as the 
arithmetic average of these quantities at neighboring integer grid 


points as 
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Vi - l (N j + VP (2.110) 

where 

N j - V (1 + ^>lj (2.111) 

When the continuity and momentum equations are solved in a coupled 
manner using an inverse procedure, four variables, F, i/>, g and 3 are 
generally treated as unknown quantities, while flow property variables 
like G, p and y are regarded as known quantities because they are 
supplied from solutions of the uncoupled equations. Since these unknown 
variables are included in the coefficients at ith level, the finite- 
difference equations. Equations (2 . 106) - (2 . 108) , are algebraically 
nonlinear. There are several ways to linearize these equations, as 
discussed in Anderson et al. (1984). Kwon and Pletcher (1981), who 
studied the effect of several linearization methods for separated flows, 
reported that among the methods considered only the Newton linearization 
with coupling of the continuity and momentum equations resulted in a 
well-behaved solution when large separation regions occurred. The 
Newton linearization with coupling is also believed to accelerate the 
convergence of the iterative solution (Keller and Cebeci, 1972; Carter, 
1978; Kwon and Pletcher, 1981). 

In the present study, the Newton linearization with coupling was 
used on the nonlinear terms in Equations (2 . 106) - (2 . 108) . The density 
and the diffusion terms were iteratively updated after solutions of the 
energy equation and the turbulence modeling equation were obtained. The 
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s c > MAX [ l« F I, 16 |/* e ] (2.113) 

j j 

-3 

and z was set equal to 1 x 10 in all the present calculations. 

The conventional way of providing the initial guess for the 
provisional value is using the value from the converged solution at the 
previous station (lagging) and 5-10 iterations were usually required for 
local convergence in compressible flows. Sometimes, 20-30 iterations 
were required near the point of separation. In the simultaneous 
viscous-inviscid interaction method, the viscous calculations take most 
of the computing time and a large number of local iterations due to the 
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linearization is undesirable. It was possible to reduce the number of 
iterations in the simultaneous viscous-inviscid interaction procedure by 
using a pseudo-time dependent approach, that is, the initial values were 
taken from solutions of the previous interaction sweep if the overall 
interaction process was reasonably converged. Convergence for the local 
Newton linearization procedure was usually obtained if the level of 
convergence error for the overall interaction process was less than 5%. 
The convergence criterion for the overall interaction procedure is based 
on the relative change of the nondimens ional mass flux m (see Equation 
(4.2)). For the first few, usually 10, global interaction sweeps, the 
initial guesses for the first provisional values were provided 
conventionally and then switched to the pseudo-time dependent approach. 
Although this remedy required additional storage for dependent 
variables, F, and G, over the interaction region, this proved to be so 
efficient that it required only 2-3 iterations at each streamwise 
station. Thus, the total computing time was reduced to only about 20% 
of that needed when lagging from the previous station is used. However, 
in the semi-inverse interaction method, the savings in computing time 
was relatively small because the inviscid calculations took the largest 
portion of the total computing time. 

Using the Newton linearization with coupling as described above. 
Equations (2 . 106) - (2 . 108) can be written as 

b j F j-l + d /j " + Vl = h / + S j 8 + C j (2.114) 

vl-i + vi + vl+i + vi - v + v + c j 


(2.115) 
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The coefficients in the above equations, which are assumed known from 
the solutions of the previous (i-1) station and previous interaction, 
are given in Appendix A for both the FLARE approximation scheme and the 
windward difference scheme. 

The appropriate boundary conditions must be given in a discretized 
form in order to complete the finite-difference system of governing 
equations. At the wall surface, Equation (2.101) becomes 

F x = = 0 (2.116) 

In order to maintain the same form as the coefficients of the blocks for 
2 < j < NJ-1, the coefficients are specified as 

A 1 = B 1 = H 1 = S 1 = C 1 = b l = S 1 = h l = C 1 = 0 (2.117-a) 

D 1 = E 1 = d l = 1.0 (2.117-b) 

At the edge of the boundary layer the continuity equation still 

holds and F^ = 1.0 should be represented in Equation (2.115) by 
choos ing 


A NJ = B NJ = E NJ = "si = S NJ = ° 


°NJ - C NJ = 10 


(2. 118-a) 
(2. 118-b) 


The procedure concerning the additional conditions required for the 
direct and inverse methods will be given in Section II. G. 3. 

b. Energy equation As was the case for the momentum equation, 

convective terms in the energy equation require appropriate treatment to 


make forward marching possible in regions of reversed flow. The FLARE 
approximation was generally used, but windward differencing was also 
used for the purpose of establishing the validity of the FLARE 
approximation. 


With the FLARE approximation, 


Cg2|F j |? c G j ■ g Vj V s ! = v N i,j + /-, G i> + c EV 2 , j+ i 


With the windward differencing, 


gW 5 G J - gS £ «.j 6,0* - «,(N 1>j+1 « n Gj) + 


( 2 . 120 ) 


where 
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N 


2,j+i 


“V * 1 - & 


+ (1 - 


_1 

Pr 




,3F | 
' 3t) [ 




(2.121-b) 


Since the continuity and momentum equations were solved first 


uncoupled from the energy equation, the velocity variables, F and ip, 
were known. Because of this, the energy equation is linear in G except 
for the density and the viscous diffusion terms which can be iteratively 
updated. However, since the diffusion coefficients in the momentum 
equation are also strong functions of temperature, the energy equation 
was solved iteratively together with the iterative solution procedure 
for the continuity and momentum equations. 

By grouping coefficients of the unknown G's, Equations (2.119) and 
(2.120) can be rewritten in tridiagonal form as 
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B.G 1 , + D.G^ + A .G 1 = C. (2.122) 

J J-1 J J J J+l J 

The coefficients are known from intermediate solutions of the continuity 
and momentum equations and the solution of the energy equation at the 
previous station and are given in Appendix B. These algebraic equations 
can be easily solved using the Thomas algorithm. 

Boundary conditions for the energy equation should be prescribed in 
a finite-difference form. When enthalpy is fixed at the wall, 


B i A i b nj a nj 0 


D i d nj c nj 1,0 


C = G 

1 W 


(2.123-a) 

(2.123-b) 


If the heat flux is given, a Taylor series expansion or polynomial 
fitting can be used to obtain the appropriate finite-difference 
expression. However, enthalpy was always specified at the wall in the 
calculations of the present study. 

c. Johnson and King turbulence model The transport equation 

for the maximum Reynolds stress given by Equation (2.75) is discretized 
with a first-order-accurate backward finite-difference approximation. 

6 e 1 = N-0 1 + N. (2.124) 

x 3 4 ' 


where 
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(2 . 125 -a) 


N. = - N 0 (1 

4 3 eq 


C, _ L 
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l-o(x)h) 


a 1 6[0.7-(y/6) m ] 


(2.125-b) 
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This equation can be rewritten simply as 


e 1 = A.e 1 ’ 1 + b. 

i i 


(2.126) 


where 


A. = 


i 1 - N^Ax 


N. Ax 
4 


B i 1 - N 3 Ax 


(2.127) 


Using an initial condition given by Equation (2.76), Equation (2.126) 
can be easily integrated along the x coordinate. During the Newton 
iteration process for the continuity and momentum equations at each 
local station, Equation (2.126) is also repeatedly solved. 


3. Method of solution 

In this section, the solution procedure for the boundary- layer 
equations is presented. As discussed in Section II. G. 2, the 
continuity and momentum equations are solved in a coupled manner and the 
energy equation and the transport equation arising in the Johnson-King 
turbulence model are solved in an uncoupled manner. The pressure 
gradient, (5, is evaluated either directly using the specified edge 
velocity (direct mode) or indirectly through the specified mass flux, m, 
(inverse mode). The parameter m is also either known explicitly from 
the interaction law (semi-inverse procedure), or determined implicitly 
during the solution procedure coupled with the inviscid scheme 
(simultaneous procedure). 
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Equations (2.114) and (2.115) can be rewritten for j = 2 to NJ-1 as 


F l-1 

D. E. 
J J 

F 1 

J 


I 

o 

[ F W 


+ 


+ 



♦1-1 

d. -1 

J J 

♦i 


0 0 

♦l+i 


H . 3+S . g+C . 
J J J 


h . (5+s . g+c . 
I J J J 


For a given £ station, these equations form a block tridiagonal system 
with each block consisting of a two by two matrix as 


[D] 1 

[A] 1 

[0] 


[Bl 2 

[D] 2 

[A] 2 

[0] 

[0] 

\ 

IB], 

x V 

A 

[A] 


\ 

N 

[0] 

\ 

V] 


J JL 



t°l I'Inj 1 d 1nj [°1nj I c >nj 


where 


A. 0 

[A], - J [»]j 

0 0 


B. 0 
J 


b. 1 
. J 


id1 j = 


D. E. 
J J 


d. -1 
. J 


H.P + S.g + C. 
J J J 


h . 3 + s.g + c . 
J J J 


|C1 : = 


The coefficients in Equation (2.130) are listed in Appendix A. 


Equation (2.129) can be solved by using an efficient form of a 
block elimination procedure which is known as the modified Thomas 
algorithm (Blottner, 1975). After elimination of the lower diagonal 
terms, [ B ] ^ , and rearrangement, the equations are reduced to a 
bidiagonal recurrence form as 


F 1 

J 


+ v + v + c ! 


i * i * * * 

= a jVi + V + s j g + c j 


(2.131) 

(2.132) 


for j - 1 to NJ- 1 . The coefficients in Equations (2.131) and (2.132) 
are given in Appendix C. 

Using the no slip boundary condition at the wall, given by Equation 
(2.117), the coefficients at j = 1 in Equations (2.131) and (2.132) 
becomes 


vc vc Vr Vc Vc Vc Vc 

A i = H i = s i = c i = a i = h i = s i = c i = 0 


(2.133) 


Then the coefficients at 2 < j < NJ-1 in Equations (2.131) and (2.132) 
can be calculated from the wall to the outer edge of the computation 
domain as discussed in Appendix C. 

With boundary conditions at j = NJ, given by Equation (2.118), 
Equation (2.132) becomes 


NJ 


= 1.0 


(2.134) 


'J ' 1 

\NJ 


^VT T ^ "F S VT ^g + C f , t 
1W iNvJ IN J 


(2.135) 


Therefore, if the pressure gradient, 3, and g(£) are known, solutions 
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for F’s and ip 1 s across the computational domain can be obtained by the 

simple back-substitution procedure starting from the edge of the 

boundary at j = NJ using Equations (2.134) and (2.135) and proceeding 

toward the wall to j = 1 using Equations (2.131) and (2.132). 

In the direct mode, the pressure gradient, f$, is determined by the 

edge Mach number, M^, which can be easily calculated from specified edge 

velocity, U^, and g(£) is defined as (2£) 2 where £ is fixed at each X 

location because £ is a function of X and U only. 

e 

In the inverse mode, if g(£) is set equal to m, the value of ip 

e 

given by Equation (2.103) can then be written in a finite-difference 
form as 

^NJ = m(I p ’ 1} (2.136) 

where 


T1 p 

I = / 6 (-)dTl 

p o p 


(2.137) 


The integral value I is iteratively updated since the energy equation 
is solved uncoupled. Also, Equation (2.135) becomes 


^NJ h NJ & + S NJ m + C NJ 


(2.138) 


Elimination of from Equations (2.136) and (2.138) results in 


B - a y m + b y 


(2.139) 




where 


a v = (I 
V P 


- 1 - 


S NJ ) 1 h NJ 


b v 


"NJ 


/ h. 


NJ 


(2.140) 


This equation describes the asymptotic relationship between the 
pressure gradient and the mass flux, m, in the finite-difference form of 
the boundary- layer equations. If the value of m is known as in the 
semi-inverse interaction method, 3 is thus determined and the solution 
can be obtained through the recurrence formulas, Equations (2.131), 

(2.132) , (2.134) and (2.135) as in the direct method. In the 
simultaneous interaction method, the value of m is regarded as an 
unknown quantity and is determined so as to satisfy an additional 
relationship between 3 and m obtained from the inviscid analysis. This 
additional relationship and the associated numerical details will be 
discussed in Chapter IV. 

Since the Newton linearization is used, the solution of the 
boundary- layer equations is repeated iteratively. Therefore, the 
coefficients in Equations (2 . 128) - (2 . 140) must be updated using the most 
recent solutions and the solutions obtained from Equations (2.131), 

(2.132) and (2.135) must converge in order to proceed to the next 
streamwise station. The convergence criterion was given by Equation 
(2.113). 

This solution procedure can be summarized as follows: 

1) Assume the initial distributions of F, if and G across the 

computational domain of the boundary layer using either lagging or 
the pseudo -time dependent approach. 
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2) Calculate the turbulent viscosity. 

3) Calculate the coefficients in Equations (2.131), (2.132) and (2.135) 
for the continuity and momentum equations. 

4) In case of the direct mode, go to step (7). Otherwise, go to step 
(5). 

5) Determine the value of m in the inverse mode (see Chapter IV). 

6) Calculate fi using Equation (2.139). 

7) Calculate the edge stream function, using Equation (2.135) in 

the direct mode and Equation (2.138) in the inverse mode, 
respectively. 

8) Calculate the solutions of F and ip across the computational domain 
of the boundary layer using Equations (2.131) and (2.132) by means 
of back-substitution. 

9) Calculate the coefficients in Equation (2.122) for the energy 
equation. 

10) Calculate the solution of G using the Thomas algorithm. 

11) Examine convergence of the solutions using Equation (2.113). If the 
solutions meet the convergence criterion, proceed to the next 
streamwise station. Otherwise, return to step (2). 
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III. INVISCID ANALYSIS 

In this chapter, the solution procedure for compressible inviscid 
flow is described. The governing equations and the appropriate boundary 
conditions are presented. The numerical methods used in solving the 
governing equations are discussed. Numerical grid generation is also 
presented. It should be noted that some of the variables used in the 
inviscid analysis may be denoted by the same symbols as used in the 
viscous analysis but the symbols may have different meanings. If this 
occurs, they will be discerned by subscripts i and v for the inviscid 
and viscous analysis, respectively. 

A. Velocity Potential Equation 

As discussed in Chapter II, viscosity effects in flows for 
sufficiently large Reynolds numbers are confined to a thin boundary 
layer near solid surfaces. As a result, the major portion of the flow 
region, which is inviscid outside the boundary layer, is governed by a 
much simpler set of equations. This reduced set of equations, called 
the Euler equations, can be obtained by neglecting viscous terms in the 
complete Navier-Stokes equations. The Euler equations consist of a set 
of first order partial differential equations, generally 4 equations 
(continuity, two momentum and an energy equations) for 2-dimensional 
flows . 

With an additional assumption of a steady, irrotational and 
isentropic flow, the Euler equations can be further simplified to a 
velocity potential equation (sometimes also referred to as the full 
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potential equation). The velocity potential formulation results in a 
single second order partial differential equation. If the flow over a 
slender body is considered, where the freestream is only slightly 
disturbed, the velocity potential equation can be simplified to form the 
transonic small disturbance (TSD) potential equation. This formulation 
provides the advantage of simplicity, especially in the application of 
the wall boundary conditions. However, the TSD equation suffers 
restrictions in applications because it is valid only for flows over a 
thin body. Recent computational evidence (Holst and Ballhaus, 1979; 
Melnik, 1981; Green and South, 1983) suggests that the velocity 
potential formulation is the most efficient one among the three (Euler, 
potential and TSD) formulations in terms of accuracy to cost ratio for a 
wide range of inviscid transonic applications. 

Isentropic flow has been assumed in developing the potential 
formulation, thus the entropy does not change across a shock wave 
according to the potential solution. The actual entropy production 
across the shock wave is proportional to the third power of the shock 
strength based on the normal component of Mach number (Liepmann and 
Roshko, 1957) as 

AS « (M 2 - 1) 3/2 (3.1) 

n 

Therefore, the assumption of no entropy change across the weak shock is 
reasonable if the normal component of Mach number, M is close to one. 
The velocity potential equation is generally regarded as a reasonable 
approximation as long as the normal component of Mach number ahead of 
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the shock wave is less than 1.3. 

The potential equation can be written and solved numerically in two 
forms known as the conservative and the nonconservative forms. South 
and Jameson (1973) developed a numerical solution scheme for the full 
potential equation in the nonconservative form using a line relaxation 
algorithm. This method was implemented into the RAXBOD code by Keller 
and South (1976) and the computational code has been used widely by many 
investigators. Newman and South (1976) reported that the 
nonconservative scheme generated a source of mass at the shock wave, 
which was caused by switching from the upwind to the central 
differencing without a proper transition operator. The strength of this 
mass source is not unique and depends on the local mesh size. A 
streamline pattern is deflected substantially by this effective mass 
source downstream of the shock wave so that the global mass balance is 
destroyed. As a result, the nonconservative scheme produces a weaker 
shock located upstream of that which would be obtained by a result of 
the Euler equations (Lock, 1981; Green and South, 1983). It should be 
noted that the fortuitous good agreement between the nonconservative 
scheme and experimental observations is due to the fact that this 
effective mass production happens to simulate the interaction of the 
shock wave and the viscous layer. 

Use of the conservative scheme maintains the conservation of mass, 
but often results in an overprediction of the shock strength (Steger and 
Baldwin, 1972; Green and South, 1983). However, there is increasing 
evidence by Melnik (1978, 1981) and Lax (1954) that a correct approach 
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should be based on the conservative scheme with viscous interaction used 
to correct the shock strength and position. In the present study, the 
velocity potential equation in the fully conservative form is used to 
describe the inviscid transonic irrotational flow. 

With the assumption of a steady isentropic inviscid flow without 
body forces or external heat transfer, the conservation laws introduced 
in Section II. A are simplified as 
continuity 


V( P V) = 0 
momentum ana energy 


( 3 . 2 ) 



JL E + V!V 

Jf-1 p 2 


( 3 . 3 ) 


According to Crocco’s theorem, a steady, inviscid, adiabatic, isentropic 
flow is also irrotational, thus permitting the use of a velocity 
potential, 0, defined as 


V = V0 


( 3 . 4 ) 


For the isentropic flow of a perfect gas, 

( 3 . 5 ) 

( 3 . 6 ) 


p 

y = constant 

o 


a = (YRT) 


i 


where a is the speed of the sound in a perfect gas. 
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1. Coordinate system 

The coordinate system chosen for the present inviscid analysis is 
presented in Figure 5. For axisymmetric external flow, an axisymmetric 
cylindrical coordinate system is used, where the abscissa, x, is 
measured along the axis of symmetry and the ordinate, r, is measured in 
the direction normal to the axis as shown in Figure 5. The velocity 
components in the x and r direction are noted as u and v, respectively. 
All the body configurations considered in this study are assumed to 
extend from the inflow to the outflow boundary. A general computational 
domain in physical coordinates looks like Figure 6. 

2. Governing equations and boundary conditions 

a. Governing equations The continuity equation can be written 
in the aforementioned coordinates as 

f^Cpur) + f^(pvr) = 0 (3.7) 

With the introduction of the velocity potential into the continuity 
equation, the potential equation is given by 

( p * x r) x + ( P^ r r) r = ° (3.8) 

b. Boundary conditions A solution of Equation (3.8) can be 
obtained by prescribing proper boundary conditions along the boundaries 
of the computational domain as shown in Figure 6. 

In the far field with a subsonic freestream, all perturbations are 
required to vanish as the flow approaches uniform conditions. 


FIGURE 5. Coordinate system for the inviscid analysis 
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FIGURE 6. Inviscid computational domain in physical coordinates 
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Therefore, the potentials at the boundaries which are placed 
sufficiently far from the source of disturbance can be assumed to be 
held at the fixed values which are determined by the uniform freestream 
condition. Assuming that the inflow boundary, noted as A-B in Figure 6, 
is placed far upstream normal to the freestream velocity, u^, the 
boundary condition is expressed in terms of the velocity potential as 

*(x Q ,r) = 0 (3.9) 

The outer boundary is also assumed to be placed sufficiently far 
from the solid wall so that its effect upon the main structure of the 
flow is negligible. Hence, the condition at the outer boundary, B-C in 
Figure 6, is also given as a constant, uniform velocity. When the B-C 
plane is parallel to the freestream, the velocity potential at r = r^ is 
given as 


0(x,r ) = u (x - x ) 
e °° o 


(3.10) 


The boundary condition at the exit of the computational domain, C-D 
in Figure 6, needs a careful treatment for convergence of the solution. 
Since the downstream boundary is positioned such that the outflow can 
always be assumed to be subsonic, a boundary condition is required. The 
downstream boundary condition at x = x^ is specified so that the 
streamwise variation of the flow is zero. 


3 2 0 

9x 2 


x=x 

e 


= 0 


(3.11) 
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The boundary condition at the inner boundary surface, D-A in Figure 
6, is the key to the viscous-inviscid interaction. Inviscid flows 
should be computed over an equivalent body, i.e., effective displacement 
surface due to the viscous effects rather than the actual body surface. 
This, however, requires repeated generation of the grid because the 
displacement thickness changes during the interaction process. In order 
to avoid this computational effort, the coupling procedure is 
implemented by using a transpiration boundary condition following the 
incompressible analysis of Lighthill (1958) who suggested that the 
effect of boundary layer displacement upon the inviscid outer flow could 
be represented by a distribution of equivalent sources on the physical 
body surface with strengths given by the streamwise growth of the 
displacement thickness. This concept has been extended to axisymmetric 
compressible flow (Gersten, 1974; Lock, 1981) as 


v 

o 


3 ~(p u r 6 ) 
ds r e e o ' 


(3.12) 


where v q is the transpiration velocity normal to the body surface, s is 
the coordinate tangential to the body surface, r Q is the radius of the 
body, and 6 is the displacement thickness. In Section II. C, the 
advantages of using the shear-layer coordinate were explained. When the 
inviscid flow is solved over the shear- layer coordinate rather than the 
physical body surface, the displacement thickness in Equation (3.12) 
should be reduced by the distance between the shear- layer coordinate and 
the actual body surface (noted as t in Figure 3). As a result, Equation 
(3.12) becomes 
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v 

o 


1 


P r 
e o 



_1 dm 

p r ds 
e o 


where 


m = p u r 6 
r e e o 


nr = p u r t 
t e e o 


(3-13) 


(3.14) 

(3.15) 


3 . Nondimens ionalizat ion 

It is convenient to use a dimensionless form of the variables in 
solving Equation (3.8). The density and velocity components are 
nondimens ionalized by the stagnation density, p Q , and by the critical 
speed of sound, a , respectively. Independent variables are normalized 
by the total length between the inflow and outflow boundary points in 
the x coordinate, L = x - x . 

3 rv 



x-x 

o 



e 




(3.16) 


m 

— m — _ t — " _ 

m “* * t 2 m +- _ * r 2 m — 

p a L t p a L 

o e o e 


m 

p a*L 2 
r o e 


where x and x are the inflow and outflow boundary position as shown in 
o e 

Figure 6. 

With the use of the above dimensionless variables, Equation (3.8) 


can be rewritten as 


(3.17) 


(p *x r) x + (p< ^r r) T = ° 

where, using Equations (3.3)-(3.6), the dimensionless density, p, can be 
expressed as 


= r i 2 


+ 0 ? 2 )] 


i/(*-D 


(3.18) 


The corresponding boundary conditions, given by Equations 
(3. 9) -(3. 13) , can be rewritten in terms of the dimensionless variables 
as 


0(0, r) = 0 

(3.19) 

0(x,r ) = u x 
e oo 

(3.20) 

a 2 *] _ o 

3x 2 _ — 0 

X=1 

(3.21) 

— _ 1 dm 

o p r ds 

(3.22) 


Hereafter, bars denoting the inviscid nondimens ional izat ion will be 
dropped, and all variables are the nondimens ional ones unless otherwise 
specified. 


4. Coordinate transformation 

The computational domain often becomes irregular in the physical 
coordinates due to the arbitrary shape of the body configuration. The 
irregular grid caused by the arbitrariness of the boundaries may 
increase the truncation error of difference schemes and sometimes 
greatly affect the stability of the solution procedure. It also takes 
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considerable effort to implement the boundary conditions on the 
irregular boundaries. The governing equations are, therefore, usually 
transformed to a suitable computational domain from the physical 
coordinates before the solution procedure is applied. Using an 
independent variable transformation, any arbitrary geometrical surface 
in the problem can be transformed to a constant coordinate line in the 
computational set so that the boundary conditions can be implemented 
with much less effort. 

A general independent variable transformation is given as 

K = £(x,r) n = n(x,r) (3.23) 

where x and r represent the physical coordinates and £ and n represent 
the computational coordinates (see Figure 7). ? = 0 and F = £ 

max 

correspond to the points x = 0 and x = 1, respectively. Likewise, n = 0 
3nd 11 max corres P ond to r = r Q and r = r , respectively. 

Applying this general transformation and maintaining the strong 
conservative form (Viviand, 1974), Equations (3.17) and (3.18) are 
transformed into the computational domain as 


+ (£ J £) , - 0 


(3.24) 


y-l 1/U-l) 

P - [! - + V V ] 


(3.25) 


where U and V are the contravar iant velocity components normal to 
constant £ and ri lines, respectively, and are expressed as 


11 / 



A? = 1 




I 


I 

I 


FIGURE 7. Inviscid computational domain in transformed coordinates 
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U = A^ + A ^ 


v = a 2 ^ + a 3 0 ti 


(3.26-a) 

(3.26-b) 


The metric quantities of the transformation are given as 


2 2 

A = 5 + 5 

1 ■’x r 


A 2 = ^x + Mr 


(3.27) 


* 2 , 2 
A 3 = ^x + \ 


j = h - K n 

x r r x 


where J is the Jacobian of the transformation. These metric quantities 
can be computed as 


J = ' Vt> 


"I 


l = J r 

T\ 


£ = - J x 

r n 


\ = ■ J r 5 "r = J X { 


(3.28) 


B. Relaxation Methods 

It is essential to know the mathematical classification of a 
governing partial differential equation in order to develop the correct 
solution procedure. Classification is easier when the potential 
equation is written in the nonconservative form as 
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(1 ' a z) *xx " ? ^xr + (1 ' ?^ rr + ^ * r = ° 0.29) 

where a is the local speed of sound. Following the standard 
classification (Anderson et al., 1984) procedure for second order 
partial differential equations, it is found that the classification 
depends upon the signs of the coefficients. This leads to the 
conclusion that the velocity potential equation becomes 

2 2 

hyperbolic where q > a (supersonic) 

2 2 

parabolic where q = a (sonic) 

2 2 

elliptic where q < a (subsonic) 

2 2 2 
where q = <b + <t> 
x r 

Transonic flows, of course, usually include all three of these flow 
categories. Flow disturbances move very quickly downstream because the 
propagation velocity is (u + a) and very slowly upstream because the 
propagation velocity is (u - a). For a disturbance in a downstream 
region to propagate upstream, it must pass around the supersonic zone. 

As a result, it is difficult to develop an efficient algorithm for 
transonic flows. Most of numerical methods for the velocity potential 
equation have been based on relaxation schemes suitable for an elliptic 
equation and adjusted for embedded hyperbolic regions because the 
supersonic region is usually relatively small in transonic flows where 
the velocity potential equation can be used. 

It should also be noted that the velocity potential equation is 
nonlinear in nature. This makes it almost impossible to use a 
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noniterative, direct solution procedure because linearization of the 
velocity potential equation will destroy the mechanism of shock 
capturing. Therefore, most of the numerical solution procedures for 
transonic flow use an iterative relaxation method with linearization. 

Relaxation methods are generally classified as either point 
iterative or block iterative. The point iterative scheme is the 
simplest and the Jacobi and Gauss-Seidel schemes are included in this 
category. These methods are very simple, but are often relatively slow 
in convergence, especially as the number of grid points increases. 
Therefore, they are seldom used in the finite-difference analysis of 
transonic flows. 

In block iterative methods, the unknown variables at a number of 
grid points are grouped together and solved simultaneously. This is 
repeated iteratively. The line-Jacobi, line-Gauss-Seidel , successive 
line overrelaxation (SLOR) and alternating direction implicit (ADI) 
methods are included in this category. These methods often provide 
faster convergence than point iterative methods. 

Iterative procedures can be accelerated by using a successive 
overrelaxation (SOR) procedure in which the values obtained from the 
standard form of any algorithm are arbitrarily modified according to the 
following format 

,n+l .n ✓ .n+1 ,rr , . 

0 = 0 + w (0 -0 ) ( 3 . 30 ) 

where w is the relaxation parameter. Here, n denotes iteration level 
and 0 n+ ^ is the most recent value of 0 calculated from the standard 


iteration algorithm, 0 n is the value from the previous iteration as 
adjusted by the previous application of this formula, and 0 n+1 is the 
adjusted (relaxed) value of 0 at n+1 iteration level. For the iterative 
procedure to converge, w must be restricted to values between 0 to 2. 
When w ranges between 1.0 and 2.0, the procedure is called 
overrelaxation. In some cases, underrelaxation is applied with w 
ranging from 0 to 1.0. 

It is known that the rate of convergence is usually sensitive to 
the choice of w and use of the optimum value of w reduces the 
computation time greatly. For a simple Laplace equation with simple 
boundary conditions, it is possible to determine the optimum u>, w 

opt 

based on the number of mesh points (Young, 1954). For complex nonlinear 
elliptic equations, it is, however, very difficult to determine the w 

opt 

in advance (Forsythe and Wasow, 1960; Ames, 1977). In such cases, 
trial-and-error numerical experimentation is often used to determine the 
value of w close to w 

opt 

The SLOR algorithm is one of the simplest block iterative 
procedures. Over the years, SLOR algorithms have proved to be reliable 
and flexible and they are still used in many numerical solution schemes 
for the potential equation (Jameson et al . , 1976; Hafez et al . , 1979; 
Chen and Caughey, 1980; Green and South, 1983). The SLOR procedure was 
used for part of the calculations in the present study and will be 
described in more detail below. 

An alternating direction implicit (ADI) method (Peaceman and 
Rachford, 1955) often provides faster convergence than the SLOR scheme. 
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In the ADI method, a partial differential operator is split into a 
sequence of an implicit operators in alternating directions, so that the 
complete iteration cycle consists of consecutive sweeps over all rows 
and columns. Since each iteration cycle consists of multiple sweeps in 
different directions, each grid point in the entire mesh can be 
influenced by every other grid point during each iteration; thus, faster 
convergence is possible. The ADI method was used in the present study 
to generate the computational grid and will be described later. The 
efficiency of the ADI method is usually dependent on flow 
characteristics and mesh sizes. However, numerous computational results 
indicate that the ADI method can be as much as 5 ~ 10 times faster than 
the SLOR method for some problems. 

The splitting of the partial differential operator into a sequence 
of one-dimensional steps as done in traditional ADI methods can also be 
formulated through an approach known as "approximate factorization". 

The terminology "approximate factorization (AF)" was used by Yanenko 
(1971) and adopted by Ballhaus and Steger (1975). In the AF procedure, 
the implicit differential operator is approximated as a product of 
factors for multi-dimensional cases. If the original differential 
operator is nonlinear, it should be linearized before factorization. 

Each factor usually requires an implicit operator in only one direction; 
thus, it involves only a simple banded matrix. Then the errors caused 
by linearization and factorization are corrected in the solution 
automatically by iterations. Some formulations of ADI schemes are often 


referred to as AF1 schemes. 


123 


A new AF scheme, AF2 , was proposed by Ballhaus and Steger (1975) in 
the study of the unsteady low frequency transonic equation in order to 
obtain faster convergence than permitted by the ADI scheme. This 
implicit AF scheme was then applied to the steady TSD equation by 
Ballhaus et al. (1978) and to the velocity potential equation by Hafez 
et al. (1979) and Holst and Ballhaus (1979). Since the AF2 scheme is 
believed to be insensitive to the moving shock instability (Holst, 

1983), this scheme is very effective for solving steady transonic flows, 
especially supercritical cases. In the subcritical case, which does not 
have supersonic regions, the ADI algorithm generally produces faster 
convergence than the AF2 scheme (Holst and Ballhaus, 1979). 

In the present study, both the SLOR and AF2 schemes were employed 
for the analysis of inviscid transonic flows. For the semi-inverse 
interaction method which employs a very simple interaction law, the AF2 
scheme was used mainly because of its high convergence speed. In the 
case of the simultaneous interaction method, the SLOR is preferred to 
the AF2 scheme since the simpler algebraic formulation of the SLOR 
reduces the human effort considerably in -manipulating the simultaneous 
interaction laws. 


C. Numerical Grid Generation 
1. Governing equations 

One of the easiest ways of establishing body-fitted coordinates is 
an analytic method which uses algebraic expressions to stretch or shear 
the coordinates. This simple procedure is useful for simple geometries 
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but is difficult to implement for complex geometries. For complex 
geometries, one of the most popular and highly developed techniques for 
grid generation establishes a new set of curvilinear coordinates from 
the numerical solution of appropriately controlled partial differential 
equations with suitable boundary conditions. 

Even though there have been a lot of different numerical grid 
generation schemes proposed (see Thompson et al., 1982), the most 
popular one has been the elliptic grid generation scheme based on the 
Poisson or Laplace equation developed by Thompson et al . (1974, 1977). 
Some of the schemes proposed more recently are based on hyperbolic 
equations (Steger and Sorenson, 1980) or parabolic equations (Nakamura, 
1982). These have the advantage that the solution can be obtained by 
simple marching from the initial plane without a time-consuming 
iteration process. However, they lack a means to control the grid 
distribution at the end of the marching plane and appear to need further 
development before they can be used in complex geometries. For the 
geometries of the present study, a scheme based on the solution of the 
Laplace equation was adequate to give smooth and nearly orthogonal 
meshes. Using a Poisson equation would have permitted additional 
control of the mesh size and skewness but this appeared unnecessary for 
the present study. 

The Laplace equation used to define the coordinate transformation 


is given by 
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5 + 5 = 0 

xx rr 


(3.31) 


ri + n =0 

xx rr 


(3.32) 


Equations (3.31) and (3.32) are transformed to the computational domain 
by interchanging the roles of the independent and dependent variables as 


A x%.v - 2B Xv + C x =0 

55 5n nn 


A r^ >- - 2B r f + C r =0 
55 5b Tin 


(3.33) 

(3.34) 


where 


A = x 2 + r 2 B = x r x + r f r C = x f 2 + r ^ (3.35) 

TIB 5 B 5 B 55 

A Dirichlet type boundary condition is prescribed at all boundaries as 

x = x(5,ti) r = r(5,n) (3.36) 


These boundary conditions are very important because they determine 
the actual grid formulation. The streamwise grid spacing along the body 
surface was chosen so that more grid points are concentrated where the 
gradients of flow properties are expected to be large. A transformation 
formula suggested by Roberts (1971) was used to produce the mesh 
clustered at some interior point, x^, along the body surface. This 
formula is listed in Appendix D. At the outer boundary, the uniform 
grid spacing was used. 

In the normal direction, grid points were concentrated near the 
body surface. This was done using an exponential stretching type 
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transformation which is also given in Appendix D. In order to start the 
relaxation scheme, initial values are needed. Initial values of x and r 
at interior points were provided by simple linear interpolation between 
boundary values. 


2. Solution procedure 

The finite-difference mesh (x(£,ti), r(£,ri)) is formed as a solution 
of Equations (3.33) and (3.34) subject to the boundary conditions given 
by Equation (3.36). For the finite-difference solution of this elliptic 
equation, a second-order-accurate central differencing scheme was used. 
Using operator notation, the finite-difference representation used for 
Equations (3.33) and (3.34) is 


L( ).. = [A. .6 r > - 2B . ,6 r + C..6 ]( ).. 

ij 1 ij IK ij U ij tui jv 'ij 


(3.37) 


where the subscript i and j represent the finite-difference mesh 
position and the finite-differencing operators are given by 

'£!< + < Ji-i.il 


8f„( ), , = 


1 


■[( ) 


( ), 


Ctr n,j 4A^An i n+ij+i v n-ij+i 


( h+ij-i + ( h-i.j-d 


(3.38) 




for the interior points. Near the boundaries, a one-sided differencing 
scheme with second-order-accuracy was used. The spacing increments A£ 
and Arj are arbitrary and were set to unity for convenience. 
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These finite-difference equations can be solved by any standard 
relaxation procedure. Although the SLOR method is usually used in grid 
generation, the ADI algorithm was used in the present study in order to 
improve the convergence rate. For a relaxation problem governed by a 
PDE in the form of L$ = 0, it is convenient to use a two- level 
correction scheme given by 

NC n + wL<£ n = 0 (3.39) 

where L is a complete differential operator, N is a linear operator 
which determines the type of iteration method, C n is the nth- iteration 
correction term which is identified as C n = - <f n , and the w is a 

relaxation parameter. L<t> n is the nth-iteration residual which indicates 
the degree of accuracy of the nth-iteration solution, tf n , to the finite- 
difference equation and is denoted by R n . Introducing an analogy where 
the number of steps is proportional to an artificial time (pseudo-time) 
coordinate t, C is considered as representing At# t . Thus, N should be 
chosen so that this process converges in time. 

In the present study, the Peaceman-Rachford ADI scheme (1955) is 
reformulated following Ballhaus et al. (1978). The N-operator is 
approximated as the product of two tridiagonal matrix factors as 

N = " a (ct " a6 ££)(“ ’ C6 Tlt ,) (3.40) 

where a is an acceleration parameter which is considered as the inverse 
of a pseudo-time step, 1/At. With the introduction of intermediate 
correction terms, f and g, Equations (3.33) and (3.34) become 
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step 1 


(a 

(a 


A 1 ? .6 rr )f n . = awLx n , 

i,j U i,j i,j 


A 1 ? , 6 rir )g n • = awLr n , 
i,j KK i,j 


(3.41-a) 

(3.41-b) 


step 2 


, r -V / n+ l n \ jr n 

(a - C. .6 )(x. . - x. .) = f . . 

i,j mr i,j i,j i,j 


, „n r w n+1 n . n 

(a - C. .6 ) (r . . - r . . ) = g. . 

i, j nn ^ i, J i, J i, J 


(3.42-a) 

(3.42-b) 


In step 1, the values of the f and g arrays are obtained by solving 
two tridiagonal matrix equations for each T| = constant line. Then the 
new values of x and y are calculated in the second step from the 
solutions of two tridiagonal matrix equations for each £ = constant 
line. These tridiagonal matrix equations can be solved easily by using 
Thomas algorithm and the actual coefficients of the resulting matrix are 
given in Appendix E. The stability analysis shows that the ADI method 
is unconditionally stable as long as 0 < w < 2 and a > 0. In the 
present calculations, convergence was obtained with a value of w equal 
to 2. 


To achieve the fastest convergence rate, it is necessary to use a 
optimum value of a which minimizes the amplification factor of the 
error. Precise estimation of the optimum a is usually extremely 
difficult. Following the suggestion of Ballhaus et al. (1978), a 
sequence of a’s is repeated with a predetermined cycle. This sequence 
is given by 


where £ is the number of elements in the sequence. The high end point 
of the sequence, a H> is effective in minimizing the high-frequency 
errors and is chosen between 0.01 and 0.05. The lower end point of the 
sequence, c* L is effective in reducing the low-frequency errors and 
recommended values are between 0 and 0.01. The values of these 
iteration parameters influence the rate of convergence and they are 
usually optimized by numerical experimentation. The typical values of 
°H used in this study are 0.02 and 0.00001, respectively, and the 

value of £ is 8. 

The converged solution is usually obtained for a 100 x 40 grid in 

50 iterations using a convergence criterion t = 1 x 10" 3 , where e is 

c c 

defined as 


z > 
c 


I Rmax‘ 


n< 


Rmax 


(3.44) 


k . 


where Rmax is the maximum residual for both x and y at the kth- 
iteration level. A typical example of numerically generated meshes is 
shown in Figure 8. 


D. Algorithms for the Velocity Potential Equation 

When the flow is subsonic, the velocity potential equation is 
elliptic. Hence a standard relaxation scheme can be applied after the 
equation is discretized with second-order-accurate central differencing 
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FIGURE 8. Grid system for the inviscid analysis 
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approximations. The von Neumann stability analysis by Holst (1983) 
shows that the central differencing scheme is stable for subsonic 
regions but unstable for supersonic regions. Therefore, modifications 
are necessary in order to provide a stable solution in supersonic 
regions . 

To deal with this problem, Murman and Cole (1971) introduced a 
type-dependent, finite-difference relaxation scheme for the transonic 
small disturbance equation. This scheme has been extended by many 
others (Steger and Lomax, 1972; Garabedian and Korn, 1972; Jameson, 

1974) to solve transonic flows using a variety of formulations. In the 
type-dependent differencing scheme, a local flow type at each grid point 
is first determined by centrally differencing the velocity potential. 

If the local flow type is subsonic, standard second-order-accurate 
central differencing formulas are used. For supersonic regions which 
are hyperbolic, first-order-accurate upwind differencing formulas are 
applied. Hence, the physical domain of dependence is correctly 
represented by the computational domain of dependence. 

For example, assuming that the main flow direction is aligned with 
the x direction, the expression for the upwind finite-difference 
representation of the second derivatives at a supersonic point (i,j) 
becomes 


* XX Ax 2 ^ 1 ^ " 20i-1 ’j + 


(3.45-a) 


(♦, : 


xr 2AxAr i , j+1 




i> j-1 


(3.45-b) 
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1 


rr 4r 2"‘i,j + l - 2 *i,j + 




(3.45 -c) 


Using a Taylor series expansion, it can be shown that an upwind 
differencing evaluation of the second derivative is first-order-accurate 
and the leading terms in the truncation error can be regarded as an 
artificial viscosity which is given by 


Ax(u 2 - a 2 )tf xxx (3.46) 

2 2 

With this artificial viscosity which remains positive when u > a , 

supersonic marching becomes stable. However, when the x-component of 

2 2 

the velocity is subsonic (u < a ) even though the flow is supersonic 
2 2 

(q > a ), the marching scheme might become unstable because the 
artificial viscosity becomes negative. This instability occurs when the 
flow is slightly misaligned with the x-axis so that the proper physical 
domain of dependence is not included in the computational domain of 
dependence. 

The rotated differencing scheme was then introduced by Jameson 
(1974) to overcome this directional difficulty. The basic idea is to 
represent the potential equation in a local stream and stream-normal 
coordinate system. Then type-dependent differencing is used along the 
local streamline coordinate while the standard central differencing is 
used along stream-normal coordinate. This produces artificial viscosity 
which is always positive where the flow is supersonic, thus eliminating 
the marching instability problem. However, the rotated differencing 
scheme has several disadvantages, too. Because it must be swept in the 
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local stream direction, difficulties might arise with general 
curvilinear meshes. Also, this scheme becomes first-order-accurate in 
both directions. 

The above methods will be referred to as artificial viscosity 
schemes since the difference molecule is adjusted in supersonic regions 
to provide the proper zone of dependence. Associated with this 
modification is a viscosity-like term in the truncation error. The same 
effect can be achieved by the artificial compressibility (or density) 
scheme. This idea has been independently developed in several different 
forms (Eberle, 1978; Harten, 1978; Hafez et al., 1979; Holst and 
Ballhaus, 1979). In the artificial compressibility scheme, an upwind 
evaluation of the density is used in supersonic regions to provide the 
upwind bias which is equivalent to what is accomplished by the 
artificial viscosity. 

In order to illustrate the idea of the artificial compressibility 
scheme which is used in the present study and to compare it with the 
artificial viscosity scheme, it is helpful to use a conservative form of 
the velocity potential equation in one-dimension, 


(p ^x ) x = ° (3.47) 

With second-order-accurate finite differencing formulas, Equation (3.47) 
is approximated as 


( p < t > ) = 6 ( p ..,6 0 .) 

xx x l+J X 1' 


(3.48) 


where 6^ and 6^ are the backward and forward difference operators for 
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the first derivative, respectively. This form is satisfactory for 
subsonic flow regions. 

As mentioned before, for the type-dependent differencing scheme, 
the difference operator is switched to upwind differencing in the 
supersonic flow region. Following Jameson (1975), this is analogous to 
the explicit addition of an artificial viscosity to the central 
differencing. The artificial viscosity is given by 

-Ax(y0 xx ) x (3.49) 

where 

V 

V = MIN [0, p(l - -£-)] (3.50) 

a 

Jameson also has shown that Equation (3.49) is equivalent to a term with 
the form of 


-Ax(vp 0 ) 
r x x x 


(3.51) 


where 


V = MAX [0, (1 - J-r)] 
9 x 


(3.52) 


When Equation (3.51) is added, Equation (3.47) can be written as 


(p Vx = 6 x (p i+iW " 6 x [v i (p i+ | ‘ p i-i } W 


(3.53) 


In subsonic regions, this scheme is second-order-accurate and centrally 
differenced. For supersonic regions, it becomes a combination of 
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second-order-accurate central differencing and first-order-accurate 
upwind differencing due to the addition of the artificial viscosity. 

In the artificial compressibility approach, the difference 
expression given by Equation (3.53) has been rearranged by Holst and 
Ballhaus (1979) to give 


(p 0 ) 

xx 


6 (p* , 6 *.) 

X 1+? X l' 


where 


(3.54) 


P i+i = (1 - V p i+± + Vi-* (3.55) 

Comparison of Equations (3.53) and (3.54) shows that the explicit 
addition of the artificial viscosity is equivalent to using a retarded 
density. The differencing becomes more strongly retarded in the upwind 
direction as the flow becomes more supersonic. If the above scheme is 
applied in both directions for two-dimensional flows, it provides the 
upwind biasing for the streamwise term in supersonic regions, thus 
giving nearly the same effect as obtained by the rotated difference 
scheme. 

One difficulty in using the artificial compressibility scheme is 
the choice of a switching function, v, because v affects significantly 
the accuracy and stability of the solution, especially, the strength of 
shock waves as supersonic regions become larger. When v is not properly 
chosen (see Equation (3.45)), a large pre-shock oscillation can be 
observed (Holst and Ballhaus, 1979), which often results in a numerical 
instability. This can be avoided by modifying the way in which v and p 
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are computed (Holst, 1983). The details will be given in a following 
section. 


E. Numerical Method 


In the present section a finite-difference representation of the 
velocity potential equation and its solution procedure are presented. 


1 . Finite-difference representation 

a. Governing equations After applying second-order-accurate 
finite-difference approximations, Equation (3.24) is written as 


6 (^) 

V j h+i.j 


V J 


= 0 


(3.56) 


where the contravariant velocity components are computed using second- 
order-accurate finite difference formulas given by 




+ 2 ( J^i+l, j 6 T|*i+l, j + 


(3.57-a) 


= Vi.J 


+ 2 ( J h,j+l 6 ^i,j+l + h,j+i 5 n^i,j 


(3.57-b) 


The values of the Jacobian, J, and the transformation metric quantities, 
A^, A^ and A^, are first evaluated at integer mesh nodes (i,j) using 


1 ^7 

•L V / 


fourth-order-accurate central finite-difference formulas. Then the 
values of the metric quantities at half integer points are obtained by 
fourth-order-accurate interpolations . 

Using the artificial compressibility scheme mentioned in the 
previous section, this formulation becomes valid for supersonic regions 
as well as subsonic regions. Instead of direct addition of the 
artificial retarded density given by Equation (3.51) which makes the 
analysis difficult due to the general £-ti coordinate system, an 
approximate implementation is achieved by using the following form 


‘ A5Cvp 5*j^$ (3.58) 

where v is still the same switching function (Holst, 1979). To simulate 
the full effect of the rotated differencing scheme, the same form of the 


artificial retarded density was added in the ri direction (Holst and 
Albert, 1979) 


i 

I 

( - Ati (vp 


Mj 

n J J 


T1 


Then the retarded density coefficients are expressed as 


- (l ' v i+M> p i +1 ,j + V.j'ira-ij 
p i,j+i - a ' + v i , j+Jt p i, j+2e-i 


(3.59) 


(3.60-a) 


(3.60-b) 


where 
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o 

II 

when 

u i+i,j 

> 

0 


( i 

when 

U i+2 , j 

< 

0 

(3.61) 

II 

o 

when 

V i, j+i 

> 

0 


1 1 

when 

V i,M 

< 

0 



The density given in Equation (3.24) can be rewritten as 

vi -Or-n 

p = U ' §^Q) (3.62) 

2 2 

where Q = A. 0*. + 2A_0 r 0 + A o 0 . In the TAIR code written by 

x 1 £ 2 fn 3 ri 3 

Dougherty et al. (1981), the density values are computed using a 
binomial series expansion of Equation (3.62). When the first four terms 
are retained, the density becomes 

p = 1 + C X Q + C 2 Q 2 + C 3 Q 3 (3.63) 

where 


C 


1 


1 

y+i 


2-Y 

2(y+i) 2 


(2-Y) (3-2Y) 
6(Y+1) 3 


(3.64) 


where X is the ratio of specific heats. 

According to Dougherty et al. (1981), this approximation is very 
accurate and saves significant computation time compared to the 
exponentiation operation. Values of the density are calculated and 
stored at the centers of the mesh cell, (i+ijj+i)? based on the 
suggestion of South and Jameson (see Holst, 1983) who found that the 
pre-shock oscillation can be reduced in this way. The values needed at 
(i + 2 >j) and (i , j+ 2 ) are obtained by simple arithmetic averaging. 
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The magnitude of the artificial viscosity, i.e., switching function 
v, is very important for the accuracy and the stability of the solution. 
Since the expression given by Equation (3.52) proved to be inadequate 
for flows with large supersonic regions, a formulation suggested by 
Holst (1979) was used in the present study: 

v if j = MAX [0, (M.,. 2 - l)C v ] (3.65) 

where M is the local Mach number and C is* a constant. With this 
1 j J V 

formulation, it is easier to select an appropriate value for v C is 

v 

usually set between 1.0 and 2.0 and should be carefully determined by 
trial -and- error depending on the strength of shock waves and relaxation 
schemes employed. The smaller values of usually produce sharper, 
less-smeared shock profiles and should be used for weak shock wave 
flows. The accuracy and stability of the solution seem more sensitive 
to the choice of in the SLOR method than in the AF2 method. This is 
believed to be due to the effect of an added artificial time-dependent 
dissipation term in the AF2 scheme. The typical value of C for a 
stable solution is 0.8 for the SLOR method and 1.2 for the AF2 scheme. 
Values of v are actually computed at the centers of the mesh cell 
( » J + 2 ) > like the density, and simple arithmetic averaging is used to 
obtain the value at (i,j). 

b. Boundary conditions The boundary conditions at the inflow 
and outer boundaries become 


0(1, j) = 0 


0(i,NJ) = « (i,NJ) 


(3.66) 
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where is the velocity potential of the undisturbed freestream and can 
be calculated easily from Equation (3.20). 

The exit boundary condition is approximated as 


£Ur _ _ r £Ur. 

1 J ^ J ; NI,j ^ J ; NI+i,j 


(3.67) 


because the £ coordinate is assumed to be almost parallel to the x 
coordinate in the vicinity of the exit plane. This approximation makes 
it much simpler to implement the boundary condition without serious 
error. The difference approximation for ^ at i = NI appearing in the 
expression of (V/J)^^ can be obtained from the boundary condition 

instead of using a one-sided differencing. From Equation (3.26-a), <p ^ 
becomes 

A, 

J - _ 

V 


U 2 

= ( r • rV 

NI, j+i 1 1 ^ NI,j+i 


(3.68) 


Consequently , 


( 


pVr % . 22 U,, 

J ) NI,j+i {pr[( J ' A x J )0 ti Aj J ]} ni . 


A 3 A 2 2 


H 


(3.69) 


Using Equations (3.67) and (3.69), the L-operator is constructed 
maintaining s econd- order- accuracy . 

A simple mathematical operation based on the coordinate 
transformation and the definition of the contravar iant velocity shows 
that the transpiration velocity boundary condition at the solid surface, 
given by Equation (3.22), can be expressed as 


v o - k 


j=i 


(3.70) 


1 A 1 


Combining Equations (3.22) and (3.70) gives 


r pVr . _ dm 

1 J J i,l d£ 


(3.71) 


This boundary condition is implemented through the control volume 
approach by assuming 


£Vr £Vr = eVr 

C J >i,lH C J i, 2 ( J } i,l 


(3.72) 


Therefore, 


r& Vr _ dm_ 

1 j J i,w z di 


- (^) 
1 J 


(3.73) 


Also, the difference approximation for at the wall appears in the 
expression for (U/J) i+ ^ ^ As for the exit boundary, Equation (3.26-b) 
is used to give 


V ^2 

= ( r • A) 

i+i , 1 3 3 K i+i,l 


Using this expression, 

A A 2 A * 

r pUr, _ A A 2 Sj . "2 dm . 

j >1+4,1 - [pr( J - a 3 j>»« a 3 d{ > I+1 


(3.74) 


(3.75) 


Equations (3.73) and (3.75) are then incorporated into Equation (3.53) 
to construct the L-operator using second-order-accurate central 
differencing formulas just as for interior points. 

“V 

As discussed before, the distribution of the mass flux, m ' , is 


assumed to be unknown in the simultaneous interaction method. 
Therefore, m must be treated as an unknown variable in developing the 


142 


general procedure for the velocity potential equation. This was 
difficult to achieve with the two step AF2 algorithm so only the SLOR 

JL 

scheme was utilized for the simultaneous interaction method and m was 
assumed to be known before each iteration for the AF2 scheme. 

2. Method of solution 

In this section, the relaxation procedures used in solving a 
finite-difference form of the velocity potential equation given by 
Equation (3.56) is presented. The present analysis employs two 
different relaxation algorithms: the SLOR and AF2 schemes. 

Just as in the numerical grid generation scheme, a two- level 
correction in pseudo-time given by Equation (3.39) was used. With this 
approach, the discrete linear operator N should be chosen so that this 
iterative relaxation process converges in pseudo-time. This often 
requires the addition of time-dependent terms to embed the steady-state 
equation in a convergent time-dependent process. With the definition of 
C n and Equation (3.30), the provisional value can be written as 

<t> n+1 = ^C n + 0 n (3.76) 

By using this expression, overrelaxation can be implemented as a part of 
the solution algorithm, not as a separate step. 

a. SLOR method The standard successive line overrelaxation 
(SLOR) procedure was applied to Equation (3.56). Even though the SLOR 
method is often much slower than the AF2 scheme, it is used for the 
simultaneous viscous - inviscid interaction method because of its 
simplicity. The simultaneous interaction method requires a relation 


between the velocity potential and the mass flux. With the SLOR scheme 
it is simple to derive such a relation because the calculation procedure 
consists of only a single sweep in each iteration. When Equation (3.57) 
is incorporated into Equation (3.56), the resulting difference molecule 
is composed of nine points as shown in Figure 9. Since the sweep 
direction for the present SLOR scheme is aligned with the main stream 
(5) direction, 0's at the downstream column, i+1, are all evaluated at 
the nth-iteration level, and 0's at the ith column are all treated as 
unknown quantities at the (n+l)th-iteration level. The 0's at the 
upstream column, i-1, are all treated as known quantities at the 
(n+l)th-iteration level. 

Using Equations (3.39) and (3.76), the resulting N is expressed as 


NC n . 
i»J 


(R. + R. .) + wR. El 1 + 6 R .6 ] C n . 
i i-i 1-15 n j n J i,j 


+ [Cross Product Terms] 


(3.77) 


The cross product terms are due to the skewness of the mesh cell and are 
represented as 




<r>i. 


4 E s 1[ V ?i ' ) i,j+i cl + 1 c ",j < 3 - 78 > 


,+l., A 2, 


where 


pA r 

R i ■ 


Rj = ( 


pA r 


J 


(3.79) 
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Odifferencing molecule for l operator 

X DIFFERENCING MOLECULE FOR N OPERATOR 


FIGURE 9. Differencing molecule for operators 
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and E. 


_ 1 


and E 


are defined as 


"s 1 ( >i,j ■ ( Vu 


,+l 


( h,i ■ < > i>J+ i 


(3.80) 


Although the correct zone of dependence is reflected on the 
computational domain through the artificial compressibility scheme, 
stability problems have been experienced near the sonic lines (where M 
is close to one) when the potential equation is solved (Jameson, 1974). 
This difficulty can usually be avoided by adding the time-dependent 
dissipation-like terms to numerical schemes (Holst, 1983). The first 
three terms in Equation (3.77) act like a time-dependent dissipation 
term (0^) in the pseudo-time analysis so that the SLOR scheme is 
usually stable for most of the transonic regime without an additional 
dissipation term. However, the SLOR scheme begins to show signs of a 
numerical instability as Mach number increases in supercritical cases. 
This instability problem at high Mach numbers might be avoided by adding 
additional time-dissipation terms (Jameson, 1974). At the present time, 
however, such an additional time -dependent dissipation term was not 
necessary when the velocity potential equation was solved coupled with 
the boundary- layer equations. 

The most important single factor determining the rate of 
convergence of the SLOR solution is the relaxation parameter, w. The 
optimum value of w is usually very close to 2, but decreases slightly as 
the Mach number increases. The typical value of w is between 1.90 and 
1.95 for moderate shock wave flows. Because of stability requirements, 
w was always set to unity in supersonic regions. 


146 


A complete ri direction operator is included in N so that this 
scheme is implicit in the B direction. On the other hand, N involves 
only the lower diagonal part of the £ direction operator, so this scheme 
is explicit in the £ direction. This implies that each grid point is 
influenced by only a single grid point to the right in the £ direction 
during one sweep so that the rate of convergence is relatively slower 
than the ADI or AF schemes. This set of N results in a tridiagonal 
matrix equation for each £ = constant line. The resulting expression 
for the tridiagonal system written for the unknown quantities at the 
ith-column becomes 


B.C n . . + D.C n . 
J i,J-l J i»J 


+ A.C n . ,, = E. 
J i.J+1 J 


(3.81) 


The coefficients of Equation (3.81), A ^ , B ^ and E ^ are functions of 
the already known C's at the adjacent columns, 0's from the previous 
iteration level, and the density. These coefficients are given in 
Appendix F. Note that only E includes dm /d£ explicitly. 

Equation (3.81) can be solved easily by eliminating either the 
upper or lower diagonal terms from the tridiagonal matrix using a 
standard Thomas algorithm. Elimination of the upper diagonal terms in 
the coefficient matrix results in the following bidiagonal recurrence 
relationships : 

C n = Pl (3.82-a) 

i,l r l 


„n , „n 

C . . = p . + q.C. . , 
1,1 J J 1 > J " J 


(3.82-b) 


The coefficients, and q ^ , are also listed in Appendix F and are 
calculated recursively from j = NJ to j = 1. After elimination of the 
upper diagonal terms, only p.^ includes dm /d£ explicitly. A closer look 

at p reveals the relation between C n and dm /d£ as 

i i,l 


„n , .dm . . .dm . . .dm’. 

C i,l P 10 P 11 d? } i-i P 12 ( d$ h P 13 d£ } i+| 


(3.83) 


If the distribution of m is known, then all C. .at the ith column can 

1 y J 

be calculated using recurrence relationships, Equations (3 . 82) - (3 . 83) . 

The SLOR solution procedure then proceeds to the next (i+1) column 
and consequently completes one streamwise sweep column by column. This 
process is repeated through multiple sweeps until convergence has been 
achieved. The evaluation of convergence is based on the degree of 
accuracy of the residual as given by Equation (3.44). The convergence 

-3 

criterion, z^> was typically set equal to 1 x 10 . A typical inviscid 

flow calculation for a supercritical transonic case requires about 
300 ~ 500 iterations. The solution procedure of the simultaneous 

interaction which treats both C ? and dm /d£ as unknown variables will 

1 y *■ 

be explained in Chapter IV. 

b. AF2 method In the approximate factorization scheme for a 
two-dimensional flow, the N-operator is usually approximated as a 
product of two factors as 


N = N 1 N 0 = L 


(3.84) 


The factors and are chosen so that their product is a close 
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approximation to L and each factor requires only simple matrix 
operations and the overall scheme is stable. Several versions of the 
AF2 scheme, which differ only slightly, (Ballhaus and Steger, 1975; 
Ballhaus et al., 1978; Holst and Ballhaus, 1979; Flores et al., 1984) 
have been suggested. The present approximate factorization scheme is 
based on the AF2 scheme which was developed by Holst (1979) and used in 
TAIR computer code. 

In the present study, the N-operator is factored as 


NC . = - -(a + 6 R.)(- a5 T 

1,J a n 2 n K 


,n 


K i K i,j 


(3.85) 


where and R^ are given by Equation (3.79). A free parameter a is 
interpreted as the inverse of an artificial time-step, 1/At, as in the 
ADI method given in Equation (3.40). This AF2 scheme can be implemented 
in a two-step format as 
step 1 


(a + 6 R.)f. . 
B J 1,J 


awL0 n . 


(3.86) 


step 2 


a6 T aXiL - 6 r R.6 r )C n . = f n . 
t> 5 K i K i,J i.J 


(3.87) 


where f. . is an intermediate result between alternative sweeps. The 
1 y J 

main difference between the present and Holst's scheme is that a regular 


right-hand coordinate system is used in the present analysis while Holst 
used a left-hand coordinate system. Therefore, some of signs in and 
in Holst's algorithm are different than those in the present 
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analysis . 

In step 1, Equation (3.86) results in a set of simple bidiagonal 
matrix equations for the f array along the T} direction. The second step 
consists of a set of tridiagonal matrix equations for the correction 
array, C, along the £ direction. The resulting sets of matrix equations 
become 
step 1 


,f n . + d.f. . 

J i.J" 1 J 



(3.88) 


step 2 


,C n , . + D.C. , 

l l-l, j i i,j 


+ A.C 
1 


n 

i+lf j 



(3.89) 


The coefficients in these equations are given in Appendix G. Equation 

(3.88) can be solved easily by back-substitution if the boundary 

condition of f at j = 1 is known. The coefficients in Equation (3.89) 

form the familiar tridiagonal matrix that can be solved by the Thomas 

algorithm after all f n, s are obtained. 

The most distinctive feature of the present AF2 scheme is that the 

difference approximation in the ti direction is split between two steps. 

This provides a 0 -type term which is thought to be helpful for the 

T) t 

convergence process of the iteration scheme as the time -dependent 
dissipation term. This also places a restriction on a sweep direction 
in both steps. The sweep direction should be in the positive ti 
direction, i.e., from j = 1 to j = NJ, for the first step and in the 

negative ti direction, i.e., from j = NJ to j = 1 , for the second step. 
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No restriction is imposed on sweep direction due to the flow direction 
in either of the two steps. 

In order to avoid the instability problem occurring near sonic 
lines, a 0^ term can be added explicitly, if necessary, in the 
appropriate factor insuring upwind differencing. In the present study, 
the 0^-type term is included by adding ± aX6^ in the second step to 
provide time-dependent dissipation in the £ direction. The difference 
direction of this term is always upwind with the flow direction as 
indicated by the double arrow notation and the sign is chosen so as to 
increase the size of the diagonal term of the matrix in Equation (3.89) 


ensuring diagonal dominance. This additional dissipation term only 
influenced stability and did not affect the values of the solution. 

The magnitude of the term is controlled by a parameter X using 

special logic suggested by Dougherty et al. (1981). The value of X is 
fixed at 0.3 in subsonic regions. In supersonic regions, X is first 
initialized depending on the shock strength. Generally, the initial 
value of X is small (~ 1.0) for small supersonic regions, moderate 
(3.0 ~ 4.0) for moderate supersonic regions, large (> 4.5) for large 
supersonic regions. Then X is updated appropriately by monitoring the 
phase of solution convergence based on the average and maximum residual, 
i.e., if the solution is converging satisfactorily, X is decreased; if 
not, X is increased. In addition to this logic, X is adjusted depending 
on how much time -dependent dissipation is required for the solution. 

This method developed by Dougherty et al. (1981) monitors the growth 
rate of the number of supersonic points, and if they grow fast, then X 


is increased; if they grow slowly, then X is decreased. 

The standard von Neumann stability analysis applied to a linearly 
simplified AF2 scheme (Holst, 1983) shows that this fully implicit 
iteration scheme is generally unconditionally stable as long as 
0 < w < 2 and a > 0. Treating the iteration level as pseudo-time, which 
can be done by considering a as 1/At, naturally suggests that fast 
convergence can be obtained with a small value of a, i.e., large time 
step. This is, however, effective for reducing only the low-frequency 
errors but not the high-frequency errors (Ballhaus et al., 1978). As 
was done in the ADI method used for grid generation, a is cycled over a 
sequence of values, given by Equation (3.43). Identifying the correct 
values for the high and low limits of a, i.e., a„ and a T , is not easy 
and is again based on trial-and-error numerical experimentation. 
Suggested values for a H and a L are 1.5 and 0.07, respectively, and the 
value of i is 8. 

Implementation of the boundary conditions outlined in the previous 
section can be done without altering the form of factors and even 
at the wall. However, a few details need attention. First, the 
boundary condition on f at j = 1 is required at the beginning of the 
first sweep. Since the intermediate variable f has little physical 
meaning related to 0 , it is not obvious how to provide a meaningful 
boundary value for f. Constructing the N-operator from the L*operator 
at the wall based on its boundary condition suggests that f^ = 0 seems 
to be a good approximation. This choice of boundary condition on f is 
also consistent with the fact that the value of f approaches zero as the 
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iteration process drives the solution to a steady state. 

Associated with the solid surface boundary condition, there is an 
additional stability condition for the present AF2 scheme. The 
parameter a has to be restricted to some minimum value in the region 
close to the solid surface when the aforementioned wall boundary 
condition is used (Holst, 1983). Usually the value of a is multiplied 
by a factor of 10 in the vicinity of the wall surface. At the exit 
boundary, the second derivative term in the £ direction in the second 
step, (6^R.6 r )C^ . 

K i K i>j 

condition . 

When the distribution of m is known, C n and 0 n+1 are determined ] 

after two sweeps of Equations (3.88) and (3.89). This process is 
repeated iteratively until convergence is achieved. The convergence 

criterion is also evaluated by Equation (3.44) and the typical value of 

-2 

was 1 x 10 . The AF2 scheme seems to converge about 10 times faster 

than the SLOR method for transonic problems. Holst and Ballhaus (1979) 
observed that even with the same level of e , the AF2 scheme provides 

i 

the solution significantly closer to the finally converged solution than 

the SLOR scheme. This is because the average residual drops more slowly I 

i 

than the maximum residual in the SLOR scheme, while both maximum and 
average residual decrease with nearly the same speed in the AF2 scheme. 

Holst and Ballhaus (1979) suggested that the reason for this | 

behavior is due to the fact that the AF2 scheme reduces all error 
components efficiently whereas the SLOR scheme treats only the high- 
frequency error components efficiently and the maximum residual is 


, is assumed to be zero according to the boundary 


highly influenced by the high-frequency errors. They also suggested 
that the root mean square error is a better criterion with which to 
compare convergence performance between different relaxation algorithms. 
Since the maximum residual is the convenient way to monitor convergence 
of the solution, the AF2 algorithm can use a larger convergence 
criterion than the SLOR scheme to provide approximately the same average 
degree of accuracy in both converged solutions. 
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IV. VISCOUS -INVISCID INTERACTION METHOD 

In the previous chapters, the solution procedures for the viscous 
and inviscid flow equations were presented. Hence, the solutions for 
each region can be obtained separately if the appropriate boundary 
conditions are provided. This zonal approach reduces the effort of 
solving the governing equations compared to the effort needed to solve 
the Navier-Strokes equations in the complete flow domain. However, in 
an interaction analysis, the solution in one region is permitted to 
influence the solution in another. In the present problem, this is 
carried out through the displacement effect of the viscous flow. This 
displacement effect causes a change in the boundary conditions for the 
inviscid flow as proposed by Lighthill (1958) and discussed in Section 
III. A. The zonal approach requires that the solutions from each zone 
match in some manner. This is usually implemented through the 
requirement that the surface pressure distribution obtained from the two 
solutions be identical. A coupling algorithm is needed to specify the 
way in which the boundary conditions for the viscous and inviscid flow 
will be altered from one iteration to the next in order to drive the 
solutions toward the matched condition. 

There have been several interaction schemes proposed and these can 
be generally classified as direct, inverse, semi-inverse and 
simultaneous (see Figure 10). In Chapter I, some of the previous 
studies employed the interaction method were reviewed briefly and were 
summarized in Table 1. If the interaction between the viscous and 
inviscid flow is weak, i.e., the viscous effect upon the pressure field 
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FIGURE 10. Coupling algorithms for the interaction method 



is small, then the classical direct interaction method can be used. As 
sketched in Figure 10-a, in the direct interaction method, the solution 
is obtained in the viscous and inviscid regions sequentially in the 
direct mode, i.e., the boundary- layer equations are solved with a 
prescribed pressure field obtained from the inviscid solution to give 
the distribution of the displacement thickness and the inviscid flow is 
solved over the prescribed displacement thickness obtained from the 
boundary- layer solution to provide the pressure field. This procedure 
is repeated iteratively until both solutions converge. This method 
usually requires underrelaxation and converges slowly. The main 
shortcoming of the direct interaction method is, however, that the 
solution of the boundary- layer equations in the direct mode leads to a 
singularity at the point of separation as discussed earlier and 
consequently destroys the interaction process. 

This difficulty can be overcome by using the inverse interaction 
method in which the role of the pressure and the displacement thickness 
is reversed from the direct method. In the inverse interaction method, 
the boundary- layer equations are solved with the prescribed displacement 
thickness and the inviscid region is solved with the prescribed pressure 
field (see Figure 10-b) . However, it is difficult to develop an inverse 
solution procedure for transonic inviscid flow. Also, it is necessary 
to use a severe underrelaxation factor even for a simple case so that 
the whole iterative procedure is very slow to converge (Carter and 
Wornom, 1975). 
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The semi-inverse interaction method proposed by Le Balleur (1978) 
offered advantages over the inverse method. In the semi-inverse 
interaction method, the boundary- layer equations are solved inversely 
and the inviscid flow is solved in the direct mode. Therefore, 
solutions including the pressure distribution in both regions can be 
obtained with the same prescribed displacement thickness. A coupling 
algorithm (denoted as interaction law in Figure 10-c) is then devised to 
update the displacement thickness in a manner which will reduce the 
difference between the pressure distributions from the two solutions at 
the next iteration. The semi-inverse strategy is shown in Figure 10-c. 
The same idea was successfully used by Carter (1979) and Kwon and 
Pletcher (1979) for incompressible interacting flows but their coupling 
algorithm was different from the one used by Le Balleur (1978). This 
method was also applied to transonic flow calculations by Whitfield et 
al. (1981), Carter (1981) and Van Dalsem and Steger (1983) and Melnik 
et al. (1983) . 

However, the coupling algorithms used in the semi- inverse methods 
to date have been rather arbitrary and lack a rigorous theoretical 
background, although Carter (1979) gave a somewhat formal justification 
for his coupling algorithm based on the von Karman momentum integral 
relation. Also, according to the von Neumann stability analysis done by 
Wigton and Holt (1981), the semi-inverse method becomes unstable for 
supersonic separated flow due to the amplification of the high frequency 
Fourier components in the solution. Such an instability has also been 
experienced in the calculation of transonic flows with large separated 
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regions in the present study. 

Since the viscous and inviscid flow regions are solved separately 
at each global iteration in the above method, only a weak coupling of 
the two regions is provided. For strongly interacting flows where the 
viscous and inviscid flow regions are locally coupled, a localized 
implicit treatment of coupling between the viscous and inviscid regions 
seems to be preferable. This type of approach which seeks a 
simultaneous solution of the the viscous and inviscid flow regions is 
classified broadly as a simultaneous interaction method (see Figure 
10-d) . This approach is further divided into the quasi -simultaneous and 
the (fully) simultaneous method depending on whether the simultaneous 
solution procedure is used partially or not. In the quasi-simultaneous 
method, the simultaneous solution procedure is used only in updating the 
boundary conditions needed for coupling, and either the boundary- layer 
or inviscid equations are solved separately. In the fully simultaneous 
interaction method, the full description of the viscous and inviscid 
flow regions are embedded together and are solved simultaneously instead 
of using a separate simple local relation. 

Most of the quasi-simultaneous and simultaneous methods developed 
to date have used integral methods for either the viscous or inviscid 
flow equations. At this time, the work of Edwards and Carter (1985) 
seems to be the only one which utilized the finite-difference method for 
both viscous and inviscid flow equations for incompressible flows. No 
solution schemes of this kind have been noted for transonic flows. In 
the present study, a simultaneous interaction method was developed for 
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the finite-difference representation of the boundary- layer and potential 
equations for transonic flow. The semi- inverse interaction method was 
also utilized in the present study. 

A. Semi-Inverse Method 

The semi-inverse interaction procedure used in the present study 
solved the boundary- layer equations in the inverse mode inside the 
interaction region and used the AF2 scheme to solve the potential 
equation for the inviscid flow. These solution procedures have been 
outlined in previous chapters. 

The coupling algorithm used in the present semi-inverse procedure 
is based on ideas introduced by Carter (1978) and Kwon and Pletcher 
(1979). In these earlier procedures, the displacement thickness 
distribution was updated based on the mismatch of edge velocities in the 
viscous and inviscid solutions employing a constant relaxation factor. 
This coupling algorithm allows use of overrelaxation which can increase 
the rate of convergence in some cases. However, it requires significant 
underrelaxation for flows with large separated regions as pointed by 
Melnik et al. (1983). This procedure is also expected to be unstable in 
supersonic separated flows according to Wigton and Holt (1981). In Le 
Balleur's semi-inverse method (1978, 1981a, 1981b), error measurement is 
based on the velocity gradient and locally optimum relaxation factors 
were applied. The preliminary calculations in the present study also 
confirmed that the use of streamwise variation in the relaxation factor 
improves the convergence rate. However, the choice of the turbulence 


1 £ 1 
±. \J ± 

model also seems to influence the convergence behavior, sometimes 
requiring careful treatment of the relaxation factor for convergence. 

In this study, a modification of the coupling algorithm was 
developed which resulted in a reduction of 10-40% in the number of 
viscous-inviscid global iterations required for a given convergence 
level. The new coupling algorithm is given as 

n+1 Ue v b 

m = m [(1 - «) + w(^) ] (4.1) 

where Ue^ is the edge velocity obtained from the boundary- layer 
solution, Uej is the tangential velocity at the surface obtained from 
the solution of the potential equation, n is the global iteration level, 
u> is a relaxation factor and b is a parameter (the modification referred 
to previously) which has been taken as 1.0 in the previous studies. The 
typical value of w used in the present calculation ranged from 0.2 to 
1.0. Letting b take on values greater than 1.0 gives greater weight to 
the local discrepancy between the viscous and inviscid solutions and has 
been found to accelerate convergence significantly without jeopardizing 
stability. In this work, values of b ranging from 1.2 to 1.7 were used. 
It should be noted that this modified form of the coupling algorithm 
generally improves the rate of convergence but not the stability 
condition of the original algorithm. 

In order to start the present semi-inverse interaction algorithm, 
an initial distribution of m is needed. This was obtained by solving 
boundary- layer equations in the direct mode using pressure data from the 
inviscid solution without the viscous effect. Usually, the adverse 
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pressure gradient from the fully converged inviscid solution is strong 
enough to cause separation in the boundary- layer calculations. Thus, 
only a few iterations were used in the inviscid solution in order to 
obtain a mild pressure gradient which did not cause boundary-layer 
separation. The distribution of dm/d£ used in evaluation of the 
transpiration velocity for the potential equation was obtained by 
applying a cubic spline interpolation method to the m distribution. 
Before applying the interpolation, m is adjusted by the distance t 
between the shear- layer coordinate and the body surface if necessary 
(see Equation (3.13)). It was also found that the computing time could 
be reduced considerably by not requiring full convergence of the 
inviscid solution at each global iteration. The number of inviscid 
solution iterations between each global viscous - inviscid iteration used 
in this work was 20 ~ 50. 

The convergence of the interaction process was determined by the 
convergence of m, i.e., 


£ > MAX [ 


n+1 n 
m - m 

n 

m 


(4.2) 


The typical convergence criterion was £ = 1 x 10 . The number of 

global iterations to obtain the above convergence varied significantly 
from case to case and will be discussed in the next chapter. 

The solution procedure of the above semi- inverse method is 


summarized as follows: 

1) Assume the position of the shear-layer coordinate, if necessary. 

2) Calculate the asymptotic solution of inviscid flow with zero 



FIGURE 11. Flow chart for the semi-inverse method 
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transpiration velocity using the AF2 scheme to obtain the initial 
guess of the pressure distribution. 

3) Using the initial pressure distribution, calculate the boundary - layer 
solution in the direct mode to obtain the initial m distribution. 

4) With the m distribution, calculate the boundary- layer solution in the 
inverse mode. 

5) If the shear-layer coordinate system was used, adjust the m 
distribution with the distance t (see Equation (3.13)). 

6) Using the adjusted m distribution, solve the potential equation for 
inviscid flow using the AF2 scheme. 

7) Examine convergence using Equation (4.2). If the solution meets the 
convergence criterion, terminate the calculation. Otherwise, proceed 
to step (8) . 

8) Update the m distribution using Equation (4.1). 

9) Return to step (4). 

The organization of this calculation procedure is also shown in 

Figure 11. 


B. Simultaneous Method 

When the simultaneous interaction method is used, the potential 
equation is solved using the SLOR scheme so that the boundary- layer 
equations and potential equation are solved along the same sweep 
direction. In the solution procedure of the boundary- layer equations, 
the local relation between the pressure gradient parameter, (S , and the 
nondimensional mass flux was given by Equation (2.139). Also, in the 
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SLOR procedure for the potential equation, the explicit local relation 

n *v 

between the correction array at the wall, C. , and dm T /d£ was 

1 3 1 1 

provided by Equation (3.83). 

By using a finite-difference approximation for dnT/d? c n can be 

i,l 

.i. 

written in terms of m . The optimum choice of the finite-difference 
expression for dm /d£ was not clear but was chosen by trial-and-error . 
The expressions used for dm /d£ are given as 


dm 

d£ 


, n+1 


.,_n+l 


. n+1 


2A£ (3m i 


" 4m i.! + m i _ 2 ) 


dm 


...n 


...n+1 


i-? 


24A£ ( '“ m i+l + 27m i 


...n+1 


,.n+l 


- 27 "i-l + m i-2 > 


(4.3) 


dm 

dt 




.t-n 


-> n 


i+i 


24AC*' - m i+2 + 27m i+l ‘ 27m i + m i-l^ 


Hereafter, A^ will be omitted because A^ was set to unity (see Section 
III. C. 2). The finite-difference expression for (d m/d£)^ is most 
important in obtaining the stable solution, because this term is 
directly associated with the transpiration velocity, v . This term was 
approximated by a second-order-accurate backward differencing formula. 
The other two terms at (i-£) and (i+i) stations usually have little 
effect on the accuracy and stability of the solution, since these terms 
appear due to the skewness of the mesh cell at the wall. They are 
approximated by fourth-order-accurate central differencing formulas. 
Especially, (dm /d£)^ +i was evaluated based on the previously known 
distribution of m at the nth-iteration level so that this term was 
treated as known quantity. 
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Since the sweep direction is along the positive £ direction, the 

values of variables for i < i-1 are already known at the (n+l)th- 

iteration level and the values of the variables for i > i+1 can be 

evaluated at the nth-iteration level. Therefore, by combining Equation 

tl _n 

(4.3) and Equation (3.83), C. can be expressed as a function of m. 

i,l l 


as 


C ”,l = Vi +X 2 


, n 


(4.4) 


The pressure gradient parameter, ft, defined in Equation (2.98) can 
be rewritten as 

dM W 


P = c 


M d£ 


(4.5) 


where 


M = u / a 
e 


(4.6) 


c m = t (1 - 1 


(4.7) 


and £^ is the transformed £ coordinate used in the analysis of the 
boundary- layer equations. Using a velocity potential, M is given as 

(4.8) 


M i + i = \ \ 


where 


Y = t^. +1 . 

i 3 i+f , 1 


i+i 

i 


(4.9) 


where J and are the Jacobian and metric quantity of the 
transformation used in the inviscid flow analysis and were given by 
Equation (3.27). 
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Applying a central difference approximation, the pressure gradient 
parameter at ith streamwise station can also expressed in terms of 
velocity potential as follows: 


L h dM* 
P I °M A$ v dCj 




^ XWl.1 


(4.10) 


M 


where M in the is approximated by values at the nth-iteration level. 
Substituting Equation (4.4) into Equation (4.10) and combining all known 
quantities, the pressure gradient parameter in the inviscid flow field 
based on the SLOR scheme can be approximated in the following form. 


,n 


n+1 n+1 


Bj. = ajtn 


+ b. 


(4.11) 


where 


a 


I 




(4.12-a) 


b I = 




i,l 




(4. 12-b) 


By nondimens ionalizing m^. in the same way as m v , and treating m fc 
in Equation (3.13) as a known quantity, the inviscid pressure gradient 
parameter, Equation (4.11), can be rewritten in terms of m as follows. 


Bj - a^m + bj 


(4.13) 
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Requiring that both pressure gradient parameters be the same, i.e., 
Sy = &i> the new value of m can be determined from Equation (2.139) and 
(4.13) as 


m 


b v - b i 
a i ' a v 


(4.14) 


When m and 8 are determined, the velocity profile in the boundary layer, 

Fj and ^ are easily calculated using Equations (2.131), (2.132), 

(2.134) and (2.135). Since the Newton 1 inear izat ion procedure was 
applied to the boundary- layer equations, the coefficients a ^ and b y must 
be iteratively updated. Therefore, m also has to be iteratively 
determined at each local streamwise station until a local boundary- layer 
solution converges. 

On the other hand, the local inviscid relation given by Equation 
(4.11) does not need to be iteratively updated because the potential 
equation analysis does not include linearization except for the density 
term which is updated at the end of each streamwise sweep. When m 

converges, C a ^ can be obtained by Equation (4.4). Before calculating 

C a 1# the relaxation factor, w, is applied to the new m. The global 
iteration process is very sensitive to the choice of w and the optimum 
value of <*> determined by trial-and-error ranged from 0.2 to 0.7. Given 

the value of C , the rest of the correction array, C n . , are 
calculated by back-substitution using the recurrence relation, Equation 
(3.82). 

The above procedure is then advanced to the next streamwise (i+1) 



station. This overall sweep procedure is repeated until the solutions 
converge. The skeleton flow chart for this simultaneous interaction 
method is shown in Figure 12. The convergence criterion used in the 
present simultaneous interaction method was also based on the 

convergence of m, which was given by Equation (4.2). The typical value 

-3 -4 

of t ranged from 1 x 10 to 5 k 10 . The initial distribution of m 

needed to start the interaction was provided in the same way as in the 
semi-inverse method. 

The calculation procedure for the simultaneous interaction method 
is summarized as follows: 

1) Assume the position of the shear-layer coordinate, if necessary. 

2) Calculate the asymptotic solution of the inviscid flow with zero 
transpiration velocity using the SLOR scheme to obtain the initial 
guess for the pressure distribution. 

3) Using the initial pressure distribution, calculate the boundary- 
layer solution in the direct mode to obtain the initial m 
distribution . 

4) With the m distribution, start the streamwise sweep of the boundary- 
layer equation and the SLOR procedure for the potential equation. 

5) Calculate the recurrence formula for C n . and the relation between 

1 > J 

and m at the ith streamwise station from the SLOR solution 
procedure . 

6) Obtain the relations between 3^ and m from the recurrence formula at 
the ith streamwise station from the boundary- layer solution 
procedure . 
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FIGURE 12. Flow chart for the simultaneous method 
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7) Calculate the new m by letting 8y = 8j. 

8) Examine convergence of m and the boundary- layer solutions. If the 
convergence criteria are met, proceed to step 9. Otherwise, return 
to step 6. 

9) Apply the underrelaxation to the converged m. 

10) Calculate the correction array, C n . with the relaxed m. 

1 > J 

11) Proceed to the next streamwise (i+1) station. 

12) At the end of each streamwise sweep, examine convergence using 
Equation (4.2). If the solutions meet the convergence criterion, 
terminate the calculation. Otherwise, return to step 4. 
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V. RESULTS AND DISCUSSION 

This chapter presents the results of the present study. The 
results are mainly divided into two categories: incompressible and 
transonic flows. First, predictions of incompressible laminar separated 
flow are presented. These were obtained as part of a preliminary study 
to test the basic performance of the present viscous-inviscid 
interaction methods. Second, the computational results for transonic 
turbulent flows with separation are presented in order to fully evaluate 
the capabilities of the present interaction schemes. 

Before the present numerical algorithm was applied to interaction 
problems, the viscous and inviscid solution procedures were tested 
separately for simple boundary layer and inviscid problems. The inverse 
solution procedure for the boundary- layer equations was verified by 
solving two laminar separated flows which have been studied numerically 
by many others; 1) linearly retarded flows studied by Howarth (1938), 
Briley (1971), Klineberg and Steger (1974), Carter (1975), Murphy 
(1977), Cebeci and Stewartson (1983), and Halim and Hafez (1984); 2) 
incompressible flows studied by Carter (1975) and Cebeci et al . (1979). 
The inverse solutions obtained by the present method for these flows 
provided good agreement with the other available solutions obtained 
based on the boundary- layer equations. The solution procedure for the 
inviscid flow was checked by comparing solutions for boattail and bump 
flows with available inviscid solutions (Wilmoth, 1977; Carter, 1981). 
The inviscid solutions obtained by the present methods (SLOR, AF2) also 
showed good agreement with the other inviscid solutions. A detailed 
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description of above comparisons is, however, not included here. 

A. Incompressible Laminar Flow 

The present viscous - inviscid interaction methods are demonstrated 
by recalculating one of flows studied by Carter and Wornom (1975). They 
analyzed a two-dimensional, incompressible, laminar, separated flow over 
a flat plate with a trough located downstream of the leading edge using 
a conventional inverse viscous-inviscid interaction method. The surface 
is prescribed as 

y = - t sech(4x - 10) (5.1) 

a 

where t is the depth of the trough as shown in Figure 13 and was set to 
a 

0.03 m in the present calculation. This provides essentially a flat 

surface far upstream and downstream of x = 2.5 m. The leading edge of 

the plate is set at x = 0 m. The Reynolds number based on freestream 

4 

conditions and the unit length is Re = 8 x 10 . In their solution 

oo 

procedure, the boundary layer was calculated in a inverse mode using 
vorticity transport and stream function equations and the inviscid flow 
was computed by the inverse Cauchy integral formulation based on small 
disturbance theory . 

This flow over a trough configuration has also been analyzed by 
many other investigators including Kwon and Pletcher (1981), Veldman 
(1981), Davis and Werle (1981), Carter and Vasta (1982b), and Edwards 
and Carter (1985). This separated flow is a good benchmark case to test 
the interaction scheme for incompressible flow, since all of the above 
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FIGURE 13. Geometric configuration of a flat plate with a trough 
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results are in good agreement. Kwon and Pletcher (1981) and Carter and 
Vasta (1982b) used the semi-inverse interaction method, and Veldman 
(1981) and Davis and Werle (1981) used the quasi-simultaneous 
interaction method. All these analyses employed the direct Cauchy 
integral formulation for the inviscid flow region. Edwards and Carter 
(1985) utilized the semi-inverse and the quasi-simultaneous interaction 
methods employing a finite-difference solution procedure for the Laplace 
equation for stream function in the inviscid flow region. All the above 
works used finite-difference solutions for the viscous region. In the 
present study, the semi-inverse and simultaneous interaction methods 
employed finite-difference methods for both the viscous and inviscid 
regions . 

The description of the present method given in Chapters II-IV was 
for the general case of compressible flow. For the incompressible 
computations, the velocity potential equation was nondimens ional ized 
using freestream conditions rather than the critical speed of sound and 
stagnation density. The rest of the solution procedure essentially 
remains the same. In the present calculation, the inflow and outflow 
boundary were set to x = -2.5 m and x = 7.5 m, respectively and the 
outer boundary was set to y = 5 m for the inviscid flow field. The 
inviscid solution was obtained with a mesh of 101 x 31 grid points in 
the streamwise and transverse directions, respectively. These mesh 
points were formed with nonuniform spacing so that a concentration of 
grid points occurred near x = 2.5 m and y = 0 m. The interaction region 
extended from x=1.0mtox=4.0m and 65 grid increments were used in 
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the streamwise direction. Inside the interaction region, the mesh for 
the boundary- layer solution used the same streamwise grid points as the 
inviscid mesh and 60 points were placed uniformly in the normal 
direction. The boundary- layer equation was started at x = 0 m from the 
uniform freestream conditions. Since the flow was assumed to be 
isothermal, the energy equation was not solved and the density and 
viscosity were assumed constant. A shear- layer coordinate was not used 
and the FLARE approximation was always used in the reversed flow 
regions . 

The same initial distribution of displacement thickness as used by 
Carter and Wornom (1975) was assumed to start the interaction procedure. 
For the semi-inverse method, the relaxation factor was set to 0.7 and 
parameter b in Equation (4.1) was set to 1.5. For the simultaneous 
method, the relaxation factor was set to 0.5. With the semi-inverse 
method, approximately 20 global iterations were required to satisfy the 

-3 

convergence criterion t = 1 x 10 (see Equation (4.2)). For the same 
convergence level, the simultaneous method converged in about 80 
iterations. However, each global iteration cycle in the semi-inverse 
method required 50 iterations for the intermediate inviscid solution to 
converge. Therefore, the simultaneous method required only about half 
the computing time required by the semi-inverse method. However, no 
effort was made to optimize the convergence process of either 
interaction scheme for this flow case. 

Some results calculated with the present interaction methods are 
compared with the predictions of Carter and Wornom (1975) in Figures 
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14-16. The results of other investigators were not plotted because all 
seem to agree very well within graphical accuracy. The predicted 
distributions of the surface pressure are shown in Figure 14. The 
present solutions are in good agreement with those of Carter and Wornom. 
Also, the semi-inverse and simultaneous methods gave almost identical 
results as expected. However, a slight difference in the prediction of 
the maximum pressure recovery at the bottom of the trough is noticeable. 
This difference is believed to be due to the different inviscid solution 
procedures used in the two interaction methods. 

Figure 15 shows the comparison of the predicted displacement 
thickness distributions. Although the present solutions agree well with 
that of Carter and Wornom (1975), the displacement thickness downstream 
of the trough is slightly overpredicted. A similar tendency was 
observed by Kwon and Pletcher (1981) and a detailed discussion about 
this can be found in their paper. The predicted skin-friction 
distribution is shown in Figure 16 along with the Blasius similarity 
solution of a flat plate. The present results show that the predicted 
values return to the Blasius solution toward the end of the interaction 
region and they are also in excellent agreement with the prediction of 
Carter and Wornom except in the vicinity of the reattachment point. 

This disagreement near the reattachment point may be due to differences 
in the grid spacing. The present mesh spacing in the streamwise 
directions was almost twice as large as that used by Carter and Wornom. 

A careful comparison of the present results with solutions of other 
viscous-inviscid interaction schemes for this laminar flow suggests that 
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FIGURE 14. Pressure coefficient distribution for an incompressible 
laminar separation bubble flow over a flat plate with a 
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trough (Re = 8 x 10 ) 
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FIGURE 15. Displacement thickness distribution for an incompressible 
laminar separation bubble flow over a flat plate with a 




FIGURE 16. Skin-friction coefficient distribution for an incompressible 
laminar separation bubble flow over a flat plate with a 
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the present interaction schemes, both the modified semi-inverse and 
simultaneous schemes, are capable of predicting the separated flows very 
accurately and they are also very consistent with each other. However, 
the present interaction methods were designed primarily for the 
transonic regime and have not been optimized for incompressible flows. 

The incompressible flow case was included to verify that the present 
algorithms do produce solutions that agree with established results in 
the incompressible limit. This was considered essential before 
proceeding to more complex transonic flow cases. 

B. Transonic Turbulent Flow 

Further evaluation of the present interaction schemes was made by 
computing transonic turbulent flows with significant separation which 
involved strong interaction between inviscid and viscous regions. 
Comparisons will be presented for two different body configurations. 

1. Boattail flow 

The first to be considered is the transonic flow over an 
axisymmetric circular-arc boattail with a solid cylindrical plume 
simulator studied experimentally by Reubush (1974). This experiment was 
conducted in the NASA Langley 16 foot transonic wind tunnel in order to 
determine the effectiveness of utilizing solid circular cylinders to 
simulate real jet exhaust plumes. Extensive measurements of surface 
pressures and boattail drag were obtained over a freestream Mach number 
range of 0.40 to 1.3 for several configurations. The boattail 
configuration is shown in Figure 17 and more detailed information can be 
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found in an article by Reubush (1974). 

In the calculations by the present methods, the computational 

domain for the inviscid flow solution, shown in Figure 18, varied with 

the body configuration; the upstream boundary was always set at 

x/d = -3.0 and the downstream boundary was set at x/d = 4.0 and 5.5 
m m 

for configurations 1 and 2, respectively. The outer boundary was always 

set at 5 d^ . The coordinate x was measured from the starting point of 

boattail and d^ is the maximum boattail diameter. The inviscid solution 

was obtained with a mesh of 101 x 31 grid points in the streamwise and 

transverse directions, respectively. These grid points were 

nonuniformly distributed so that about a third of total gird points were 

placed along the boattail region and a concentration of mesh points was 

placed at the boattail-sting juncture. The interaction region was 

assumed to extend from x/d = -1.0 to x/d = 2.5 for the configuration 1 

mm 

boattail and from x/d = -1.0 to x/d = 3.5 for the configuration 2 

mm 

boattail. Inside the interaction region, the same streamwise grid 
spacing was used at the body surface for both the inviscid and viscous 
solutions and 61 ~ 65 grid points were used. In the transverse 
direction, 70 boundary- layer grid points were nonuniformly distributed 
based on Equation (2.105). The ratio between two adjacent grid spacing, 
K, was 1.05. 

A velocity profile at the upstream boundary was obtained from the 
direct boundary- layer solution starting from the leading edge 
(x/d = -6.64) using the pressure obtained from the undisturbed inviscid 
solution. This inflow velocity profile was not updated during the 
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FIGURE 18. Computational domain for the inviscid solution of the 
boattail flow 
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interaction process because the pressure distribution upstream of the 
interaction region obtained from the interacted inviscid solution 
remained nearly unchanged. The effect of the cone-cylinder forebody was 
neglected in the present calculations. The initial distribution of the 
mass flux, m, was obtained by the boundary- layer solution in the direct 
mode using the pressure data provided by the undisturbed inviscid 
solution of the AF2 scheme with only 10 iterations. The distribution of 
m downstream of the interaction region was iteratively calculated by 
assuming that dm/dx in that region decreased linearly from the value at 
the interaction boundary to zero at the exit boundary. 

In every turbulent flow calculation for this configuration, 

transition from laminar to turbulent flow was assumed to occur near the 

leading edge at x/d^ = -6.0 and the algebraic Cebeci-Smith and Johnson- 

King turbulence models were used. When the Johnson-King model was used, 

it was initiated at x/d = -1.1 after starting with the Cebeci-Smith 

m 

model. In the figures to follow, the results of the simultaneous and 
semi-inverse methods obtained with the Johnson-King turbulence model 
will be shown with the results of the simultaneous method obtained with 
the Cebeci-Smith model. Since the results of the semi-inverse 
interaction methods obtained with the Cebeci-Smith model were always in 
good agreement with those of the simultaneous method obtained with the 
same model, they will not be displayed. In the region of reversed flow, 
the FLARE approximation was always used. 

Even though Reubush (1974) performed extensive measurements for a 
wide range of Mach numbers, Reynolds numbers and boattail 
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configurations, only a few representative cases were calculated in the 
present study. The first comparison was made for the configuration 2 
boattail case at M =0.8 and Re = 1.22 x 10^, which resulted in a 

oo oo 

fully attached flow over a smooth boattail surface without the 
appearance of a shock wave. Because of the smooth body surface, the 
shear-layer coordinate was not used in this case. With the semi-inverse 
method, the relaxation factor was set to 1.2 and the value for parameter 
b in Equation (4.1) was also set to 1.5. The solution converged nicely 

-3 

in 30 iterations to E = 1 x 10 . Each global iteration cycle in the 

semi-inverse method required about 90 seconds on a Perkin-Elmer 3240 
computer, which is similar to a VAX-11/780. For the simultaneous 
interaction method, a relaxation factor of 0.75 was used and 110 
iterations were required for the same convergence criterion. 
Approximately 20 seconds were required to complete one iteration cycle. 
As a result, both interaction methods took about the same overall 
computing time for this case. 

The predicted distribution of surface pressure is shown in Figure 
19 and is compared with the solutions of Wilmoth (1977) and Swanson et 
al. (1983). As discussed earlier, the solutions of Wilmoth (1977) were 
obtained using the direct viscous -inviscid interaction method with 
experimentally determined separation and reattachment points. Swanson 
et al. (1983) solved the Navier-Stokes equations with the composite 
velocity formulation and used an algebraic Cebeci-Smith turbulence model 
with a relaxation formula of Shang et al. (1976). The comparison shows 
that the present solutions with the Johnson-King model are in the best 
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FIGURE 20. Skin-friction coefficient distribution for a transonic 
turbulent flow over a boattail (Config. 3, = 0.8, 
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agreement with the experimental data. The results of the Cebeci-Smith 
model in the present analysis and Navier -Stokes solution by Swanson 
et al. (1983) show a slight overprediction of the pressure in the 
vicinity of x/d = 1.7. The predicted skin friction for this flow is 
given in Figure 20. Unfortunately, experimental data were not available 
for comparison. As in the comparison of the pressure distribution, the 
predictions from both the simultaneous and semi-inverse method with the 
Johnson-King model are indistinguishable. However, the predictions of 
the Cebeci-Smith model differ from the others upstream and again 
downstream of the boattail. The Johnson-King model predicts the lower 
skin friction upstream of the boattail but higher skin friction in the 
downstream region. The same tendency was observed in every other 
calculation in the present study. This difference is nearly insensitive 
to the choice of the starting point for the Johnson-King model. 

The next comparison is for the separated flow over the 

configuration 1 boattail, the most steep boattail tested by Reubush 

(1974), at the subcritical conditions of M =0.7 and Re = 1.16 x 10^. 

Experimental data indicated that a pressure plateau was formed over the 

last 30% of boattail followed by a trailing edge compression without the 

appearance of a shock wave. The separation point determined from oil 

flow visualization studies by Abeyounis (1977) on the same flow model 

was x/d = 0.51. 
m 

In the present calculations, the relaxation factor was set to 0.35 
and 0.50 for the semi-inverse and simultaneous interaction methods, 
respectively. The convergence behavior of the iteration process 
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depended on the turbulence model used. In the semi-inverse interaction 
method with the same convergence criterion as used in the previous case, 
the Cebeci-Smith model required about 70 ~ 80 iterations and the 
Johnson-King model required 100 ~ 120 iterations. In the simultaneous 
interaction method, the turbulence models used had little effect on the 
convergence behavior and the solutions converged in 170 ~ 200 
iterations. In terms of total computing time, the simultaneous method 
required only about 40 ~ 50% of that needed in the semi-inverse method. 

Figure 21 compares the surface pressure distribution predicted by 
the present interaction schemes with experimental data. Also shown in 
Figure 21 are the Navier-Stokes solutions obtained by Swanson et al . 
(1983) and the viscous - inviscid interaction results of Wilmoth (1977). 
Figure 22 shows the calculated variation of the skin-friction 
coefficient and the Navier-Stokes solution of Swanson et al. (1983). No 
experimental skin-f r iction data are available at the present time for 
comparison. The present solutions obtained with the algebraic Cebeci- 
Smith model significantly overpredict the surface pressure variation in 
the vicinity of the boattail-sting juncture. The predictions of the 
separation and reattachment points (x/d^ = 0.64 and 0.88, respectively) 
by the Cebeci-Smith model are also very poor. The Navier-Stokes 
solutions obtained with the Shang-Hankey turbulence model provide better 
predictions than the present solution obtained with the Cebeci-Smith 
model. The prediction of surface pressure in the reversed flow region 
is, however, still much higher than the experimental data, and the 
solution does not seem to capture the pressure plateau which is 
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characteristic of the extensive separation. Also, the prediction of the 
separated region is not much better than present prediction obtained 
with the algebraic Cebeci-Smith model. The interaction results by 
Wilmoth show better agreement with the experimental data than the above 
results and clearly exhibit the general characteristics of the pressure 
plateau but show a rapid pressure increase in the second recovery 
region. 

On the other hand, the present prediction based on the Johnson-King 
turbulence model generally agrees very well with the experimentally 
observed data. Especially, the prediction of the separation point 
(x/d^ = 0.56) is much better and the predicted reversed flow region is 
almost twice as large as given by the other solutions. Since the 
experimentally observed data for the reattachment point are not 
available, it is difficult to argue which prediction is better. 

However, based on the view of the other experimental data considered in 
the present study, the reattachment point usually occurs near the point 
of maximum pressure. On this basis, the prediction of the reattachment 
point by the Johnson-King model is believed to be better than the other 
solutions . 

Figures 23-25 show the distributions of the displacement thickness, 
maximum Reynolds stress and the mean velocity profiles calculated by the 
Cebeci-Smith and Johnson-King models for the same separated boattail 
flow. Experimental and other numerical data were not available for 
comparison. The displacement thickness predicted by the Johnson-King 
model is larger by as much as 30% than that predicted by the Cebeci- 
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FIGURE 23. Displacement thickness distribution for a transonic 
turbulent flow over a boattail (Config. 1, M =0.7 




FIGURE 24. Maximum Reynolds shear stress distribution for a transonic 
turbulent flow over a boattail (Config. 1, M = 0.7, 
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Smith model. This is consistent with the smaller pressure recovery 
predicted by the Johnson-King model downstream of the boattail. 

The most striking difference between the two models can be found in 
the comparison of the calculated variations of the maximum Reynolds 
stresses displayed in Figure 24. The Cebeci- Smith model predicts a 
sharp peak in -u'v' , which suggests that the turbulence field is 
strongly dependent upon the local flow conditions. On the other hand, 

-u v predicted by the Johnson-King model shows more gradual increase 
over a longer streamwise distance, which is responsible for the thick 
displacement thickness and smaller pressure recovery. From the velocity 
profiles shown in Figure 25, a major difference can be observed in the 
inner part of the boundary layer: the Johnson-King model predicts much 
larger separation bubble than the Cebeci-Smith model. 

The results for the boattail case presented in Figures 19-25 were 
obtained using a body surface oriented coordinate system for the viscous 
calculation. Guided by these results, the calculations were repeated 
using two different shear-layer coordinate systems. First, a dividing 
streamline O = 0) from the previous results was chosen as the shear- 
layer coordinate. A second choice of shear- layer coordinate was made so 
as to align the coordinate with the displacement surface in the reversed 
flow region. The coordinate then asymptotically approached the body 
surface. A comparison of the calculations using both shear-layer 
coordinates with the previous results obtained using the body-oriented 
coordinate shows differences too small to be drawn on the figures, and 
suggests that the use of the shear- layer coordinate system provides no 
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important advantage in the calculation. 

From viewing the results of this subcritical flow case, it becomes 
obvious that the new simultaneous viscous-inviscid interaction method 
can usually predict flows with pressure driven separation more 
economically than the semi-inverse method, for the same degree of 
accuracy. It was also found that the accuracy of the solutions is more 
strongly dependent on the turbulence modeling than the interaction 
algorithm. The Johnson-King turbulence model produced much better 
predictions than the algebraic Cebeci-Smith model for separated flows 
with viscous-inviscid interaction. 

2. Bump flow 

The next case to be presented is the shock interacted flow field 
about an axisymmetric circular-arc bump studied experimentally by 
Bachalo and Johnson (1979). This flow served as the primary test case 
for the present study because of the availability of extensive 
experimental data. This flow configuration, shown in Figure 26, 
consists of an axisymmetric circular-arc bump attached to a circular 
cylinder placed parallel to the flow direction. The cylinder, 15.2 cm 
in outside diameter, extends 61 cm upstream of the bump leading edge. 
The bump has a chord of 20.32 cm and is 1 . 9 cm in thickness. The bump 
leading edge is smoothed by a small fillet, of radius 18.3 cm, to 
prevent separation at the leading edge. The model was initially tested 
at a freestream Mach number, M = 0.875 and freestream Reynolds number, 
Re = 1.36 x 10 7 in the NASA Ames 2x2 foot transonic wind tunnel. 

oo 


BUMP GEOMETRIC PARAMETERS 
= 61.00 cm 

C = 20.32 cm 

d = 15.34 cm 

Rf = 18.30 cm 

Rg = 28.05 cm 


FIGURE 26. Geometric configuration of a bump 
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Later, Horstman and Johnson (1984) reported the measurements for 
the same body configuration obtained in the 6 x 6 foot supersonic wind 
tunnel over a range of freestream Mach numbers from 0.4 to 0.94. This 
experiment was carried out in order to evaluate the importance of the 
interference effect of the wind tunnel walls. For M = 0.875 both 
results agreed well on general aspects of the flow including the 
distribution of the surface pressure. However, a slight difference was 
noticeable in the location of the shock wave; the shock wave in the 
larger wind tunnel was located about 0.01 chord length upstream of that 
measured in the smaller wind tunnel. Also, the pressure recovery 
downstream of the shock in the larger wind tunnel was slightly smaller 
than that observed in the smaller wind tunnel as shown in Figure 27. 

Their experimental data can be summarized as follows. For a Mach 
number up to 0.8, the flow was subcritical (no shock wave appears) and a 
small separation bubble was present near the bump-sting juncture. The 
locations of the separation and reattachment points were somewhat 
insensitive to the Mach number in subcritical cases. As the Mach number 
passes 0.8, the flow becomes supercritical and the shock wave forms 
approximately at x/C =0.65 (x is referenced to the forward intersection 
point of the arc of the bump with the cylinder and measured parallel to 
the cylinder and C is the chord length of the bump) . As the Mach number 
is increased further, the location of the shock wave moves very slightly 
downstream and a large pressure plateau is formed in the vicinity of the 
bump-sting juncture followed by the smooth compression after the 
juncture. The separation bubble also begins to grow very rapidly; the 
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separation point moves upstream toward the shock wave location and the 
reattachment point moves farther downstream. For example, at 

= 0.875, the shock wave is located approximately at x/C =0.67 and 
the separation and reattachment points are located at x/C - 0.69 and 
x/C - 1.17. The extent of the separated region was determined by oil- 
flow visualization. Considerable uncertainty was reported to exist in 
interpretation of these observations, especially the reattachment point. 

Johnson et al. (1982) performed numerical computations using the 

Navier-Stokes equations with the algebraic Cebeci-Smith and Wilcox- 
2 

Rubesin (k-w ) two-equation turbulence models. Their results indicated 
that the shock location and the separation point were predicted 
substantially downstream of the experimentally observed position and the 
more sophisticated two-equation model did not provide significantly 
better predictions than the simple algebraic model. Horstman and 
Johnson (1984) used the k-£ turbulence model developed by Jones and 
Launder (1972) in their Navier-Stokes solutions and obtained a slight 
improvement in the predictions of the shock location and the separation 
point over the results of Johnson et al. (1982). Most of the 
improvement was, however, due to the new outer boundary condition: they 
used the freestream outer boundary condition instead of the solid wall 
boundary condition based on the experimental data which indicated that 
the flow was almost free of the wind tunnel wall effects. However, this 
resulted in an overprediction of surface pressure downstream of the bump 
trailing edge region. Recently, Johnson (1985) solved the same bump 
flow using the Navier-Stokes equations with the Johnson-King turbulence 
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Comparison of measured pressure coefficients with 
predictions for inviscid flow for a transonic turbulent flow 
over a bump (M^ = 0 . 875 ) 
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model. The overall predictions were in much better agreement with the 
measurements than the above Navier-Stokes solutions based on the 
algebraic and two-equation models. 

The semi-inverse viscous - inviscid interaction approach was applied 
to the same flow by Carter (1981) and Carter and Vasta (1982a). The 
solution was obtained with the algebraic Cebeci-Smith model and predicts 
the shock location better than the previous Navier-Stokes solutions but 
the pressure plateau was not captured at all and the pressure at the 
bump trailing edge region was substantially overpredicted. It was also 
reported that the predictions were improved by reducing the Clauser 
constant of the Cebeci-Smith model by half (Carter, 1981) and by using 
the relaxation formula of Shang and Hankey (1975) (Carter and Vasta, 
1982a). 

In the present calculations, 101 x 41 mesh points were nonuniformly 
distributed over the inviscid computational domain which extended from 
x/C = -2.0 to x/C =3.0 in the streamwise direction and 2.5 chord 
lengths in the transverse direction. Preliminary calculations were 
performed with several different computational domains and mesh sizes. 
The computational domain and the number of grid points were reduced to 
the smallest values that still provided mesh independent interaction 
solutions. About 40% of the mesh points along the body surface were 
placed over the bump (0 < x/C < 1.0) with the greatest concentration 
near x/C = 0.65. The smallest mesh spacing was Ax/C = 0.015. The 
viscous-inviscid interaction region extended from x/C = -0.5 to 
x/C = 1.5 and occupied 61 mesh points in the streamwise direction. For 
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the viscous mesh inside the interaction region, the inviscid mesh 
spacing along the wall was used in the streamwise direction and 70 
points were distributed nonuniformly across the boundary layer with 
K = 1.09 (see Equation (2.105)). 

The incoming velocity profile for the boundary- layer analysis was 
again obtained from the direct solution of the boundary- layer equations 
starting from the nose of cylinder (x/C = -3.0) using the pressure 
obtained from the undisturbed inviscid solution. The initial 
distribution of m was provided by the direct solution of the boundary- 
layer equations with the pressure data obtained by the undisturbed 
solution of the potential equation using the AF2 scheme with 10 
iterations. The distribution of m downstream of the interaction region 
was obtained in the same way as in the boattail flow calculations. 
Transition from laminar to turbulent flow was made at x/C = -2.8. 
Calculations were made using both the Cebeci-Smith and Johnson-King 
turbulence models. The switching point for the Johnson-King model was 
x/C = -0.6. 

The convergence behavior of the viscous -inviscid iterative 
procedure was found to be very sensitive to the choice of turbulence 
model and the freestream Mach number. This implies that convergence of 
the interaction process is strongly dependent on the size of the 
separated flow region. The relaxation parameter was decided by a trial- 
and-error method. In subcritical cases which resulted in small 
separation regions, a relatively large value of the relaxation factor 
was used regardless of the turbulence models (0.7 ~ 1.0 for the semi- 
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inverse method and 0.6 ~ 0.75 for the simultaneous method) and the 

solutions usually converged very well for both interaction methods 

(35 ~ 50 iterations for the semi-inverse method and about 100 iterations 

for the simultaneous method). When the semi-inverse method was used for 

supercritical cases with the Mach number less than 0.9, a relaxation 

factor of 0.4 ~ 0.6 was used with the Cebeci-Smith model and convergence 
-3 

to e = 1 x 10 was achieved in 100 ~ 150 iterations. Each global 
iteration cycle in the semi-inverse method required about 2 minutes. 

When the Johnson-King model was used in the semi-inverse method, a 
smaller value of the relaxation factor (0.2 ~ 0.3) was necessary to 

obtain convergence. However, the solutions converged to only about 

-2 

e = 1 x 10 after 40 ~ 50 iterations. It has not been possible to 

reduce the residual further due to the occurrence of oscillatory 

behavior in the convergence pattern. For M =0.9, the solution by the 

semi-inverse method diverged even with a small relaxation factor (0.1). 

On the other hand, with the simultaneous method for the supercritical 

cases with the Mach number up to 0.9, the convergence was still obtained 

with a relatively large value of relaxation factor (0.4 ~ 0.7) 

regardless of the turbulence model and 150 ~ 250 iterations were needed 
-3 

for e = 1 x 10 . A typical calculation with the simultaneous 

interaction method requires 60 ~ 80 minutes on a Perkin-Elmer 3240 
minicomputer . 

Before the interaction method was applied, the present numerical 
algorithm was carefully tested especially for the inviscid flow 
calculation. In order to obtain a meaningful comparison of the semi- 
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inverse and simultaneous interaction methods, the AF2 and SLOR schemes 
must provide nearly the same inviscid flow solutions. For subcritical 
cases, the undisturbed inviscid solutions of the SLOR and AF2 schemes 
were found to be almost identical. When the artificial compressibility 
approach is used in the solution procedure of the full potential 
equation for supercritical cases, it is already known that the 
prediction of the shock strength is very sensitive to the magnitude of 
density biasing. Therefore, the density biasing magnitude must be 
decided carefully. This was done by numerical experiments and 
comparison with the results of the other inviscid code and available 
solutions in the literature. Figure 27 compares the undisturbed 
solutions of the SLOR and AF2 schemes for = 0.875. Also shown in 

Figure 27 is the inviscid solution of Carter (1981) for comparison. 

Both solutions are in good agreement except that the SLOR solution shows 
a slightly stronger shock than the AF2 solution. This difference 
disappears when the viscous effect is included. 

Figures 28-31 compare the calculated surface pressure distributions 
with the experimental data taken from the 6x6 foot wind tunnel at 
= °-6, °-8, 0.875 and 0.9 by Horstman and Johnson (1984). The 
Navier-Stokes solutions obtained by Johnson (1985) using the Johnson- 
King turbulence model are also compared for all the above Mach numbers. 

^oo ~ 0.875, the experimental data of the 2x2 foot wind tunnel 
(Bachalo and Johnson, 1979) and the Navier-Stokes solutions obtained by 
Horstman and Johnson (1984) based on the Cebeci-Smith and k-e turbulence 
models are also displayed. The Cebeci-Smith model used in the present 
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interaction scheme overpredicts the pressure in the separated region and 
does not capture the pressure plateau characteristics at all. The same 
trend is also observed in the Navier-Stokes solutions of Horstman and 
Johnson (1984) using the same Cebeci-Smith turbulence model, but the 
present predictions of the shock position seem to agree slightly better 
with the experimental data than the Navier-Stokes solutions based on the 
same Cebeci-Smith turbulence model. The overall predictions of the 
present interaction method using the Johnson-King model are observed to 
agree very well with the experimental data and the Navier-Stokes 
solutions of Johnson (1985) with the same Johnson-King turbulence model. 

The present semi-inverse and simultaneous interaction methods again 
give almost identical final solutions except for M =0.9 for which the 

oo 

semi-inverse method failed to converge. The shock position is 
reasonably well predicted and the pressure plateau behavior is well 
described by the Johnson-King turbulence model for all Mach numbers. 
However, in the supercritical cases, the shock strength is slightly 
underpredicted and the shock location is predicted upstream of the 
measured data by both the semi-inverse and simultaneous interaction 
method. The difference between the solutions obtained by the present 
interaction method and the Navier-Stokes solutions based on the same 
Johnson-King model becomes larger as the Mach number is increased. One 
interesting observation can be made for the M = 0.875 case. While the 

oo 

Navier-Stokes solutions of Johnson (1985) agreed very well with the 
experimental data of the 2x2 foot wind tunnel, the present predictions 
agree better with the 6x6 foot wind tunnel data. Considering that the 
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FIGURE 29. Pressure coefficient distribution for a transonic turbulent 
flow over a bump (M =0.8) 
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FIGURE 31. Pressure coefficient distribution for a transonic turbulent 
flow over a bump (M =0.9) 
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wind tunnel wall effect should be much smaller in the larger wind 
tunnel, it could be argued that the present interaction solutions 
provide the best agreement with the experimental data. 

In the present calculations for this flow case, solutions were 
obtained using both the FLARE approximation and the windward 
differencing scheme in the boundary-layer equations. Figure 31 compares 
the pressure coefficient predicted using the FLARE approximation with 
the predictions using the windward differencing for the flow at 
= 0-9. This flow has the largest separation bubble among flows 
considered in the present study. The two solutions gave nearly 
identical results. The convergence behavior seems to be insensitive to 
the use of the FLARE approximation and the windward differencing scheme. 

The skin-friction coefficient predicted by the present interaction 
scheme using the Cebeci-Smith and Johnson-King models for M =0.6, 

0.875 and 0.9 are shown in Figures 32-34. Experimental data and other 
numerical solutions for skin friction are not available at the present 
time. Again, there is little difference between the solutions of the 
semi-inverse and simultaneous method. For M = 0.6, the difference 

oo 

between the predictions of the Cebeci-Smith and Johnson-King turbulence 
models is very small and the predicted separation and reattachment 
points are in good agreement with the experimental data. 

For = 0.875, the difference between the Cebeci-Smith and 
Johnson-King turbulence model becomes evident. The Cebeci-Smith model 
results in a much smaller separation region and the separation and 
reattachment points predicted by the Johnson-King model are seen to be 
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FIGURE 32. Skin-frict ion coefficient distribution for a transonic 
turbulent flow over a bump (M =0.6) 
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FIGURE 33. Skin-friction coefficient distribution for a transonic 
turbulent flow over a bump (M = 0.875) 
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FIGURE 34. Skin-friction coefficient distribution for a transonic 
turbulent flow over a bump (M =0.9) 
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FIGURE 35. Comparison of separation and reattachment locations for a 
transonic turbulent flow over a bump 
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in fairly good agreement with the measurements. As was observed in the 
calculation of the boattail flow, the skin friction predicted by the 
Johnson-King model shows a more rapid increase after the reattachment 
point than the prediction by the Cebeci-Smith model. In Figure 34, the 
predictions of skin friction by the FLARE approximation and windward 
differencing are compared for = 0.9. Some differences can be 
observed in the region before the reattachment point but overall, the 
agreement is reasonably good and the predictions of the separation and 
reattachment points agree very well. Therefore, it was felt that the 
use of the FLARE approximation was justified for the flows considered in 
the present study. 

Figure 35 compares the predicted separation and reattachment points 
with the experimentally observed values and values from the Navier- 
Stokes solutions obtained using the Jones -Launder (k-e) turbulence model 
(Horstman and Johnson, 1984) and Johnson-King turbulence model (Johnson, 
1985). The predictions obtained by the present simultaneous interaction 
method with the Johnson-King model is generally in good agreement with 
the measurements. The reattachment point is predicted slightly 
downstream of the measured value. The predictions by the Cebeci-Smith 
model are relatively good in lower Mach numbers, but the disagreement 
with the measurements becomes larger as the Mach number is increased 
above 0.85. On the other hand, the Navier-Stokes solutions by Johnson 
(1985) with the Johnson-King model are observed to significantly 
underpredict the separation point compared to the measurements and the 
present results with the same model. It is also interesting to note 
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that the Navier-Stokes solutions with the k-e two-equation model failed 
to prfedict separation when the freestream Mach number was less than 0.8. 

Note that the experimental data of the M = 0.875 case shown in the 

oo 

figures to follow are the measurements taken in the 2x2 foot transonic 
wind tunnel. Measurements from the 6x6 foot wind tunnel are not 
available except for the surface pressure distribution at the present 
time. The predicted and measured displacement thickness distributions 
of the boundary layer for M~ = 0.6 and 0.875 are compared in Figures 36 
and 37. None of the predictions agree really well with the measurements 
over the full extent of the flow. Discrepancies are especially evident 
after the bump-sting juncture (x/C = 1.0). The difference becomes 
larger at the supercritical Mach numbers. The present solutions based 
on the Johnson-King model predict the largest displacement thickness. 
This is also consistent with the fact that the pressure recovery along 
the bump trailing edge predicted by the Johnson-King model in the 
present study is relatively the smallest. 

Figures 38 and 39 show the comparisons of the predicted maximum 
Reynolds stresses and the experimental data for M =0.6 and 0.875. In 

oo 

the comparison of the Reynolds shear stress of the boundary- layer 
calculations with the measurements, the effect for the rotation of 
coordinate along the body surface may be significant and must be 
corrected. The experimental data shown in Figures 38 and 36 are the 
values recalculated by Johnson (1985) from the measurements of Bachalo 
and Johnson (1979) by performing a coordinate transformation from the 
rectangular measurement coordinates to shear-layer coordinates. This 
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FIGURE 36. Displacement thickness distribution for a transonic 
turbulent flow over a bump (M =0.6) 
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FIGURE 37. Displacement thickness distribution for a transonic 
turbulent flow over a bump (M = 0.875) 
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FIGURE 39. Maximum Reynolds shear stress distribution for a transonic 
turbulent flow over a bump (M =0.875) 
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comparison shows the most striking difference between the algebraic 
Cebeci-Smith and Johnson-King model. The Cebeci-Smith model shows more 
rapid increase and abrupt decrease in -u’v’^, while the Johnson-King 
model provides more accurate prediction near the shock wave and smoother 
development across the flow field. One thing to notice is that the 
present predictions with the Johnson-King model indicates a slower decay 
downstream of the reattachment point than do the Navier -Stokes solutions 
obtained by Johnson using the same model. This slower decay is believed 
to be the cause of the rapid increase in the skin-friction coefficient 
after the reattachment point predicted by the Johnson-King model in the 
present method. In his Navier-Stokes solutions, Johnson (1985) found 
that this slow decay in the maximum Reynolds stresses was due to the 
fact that the value of a(x) in Equation (2.73) grew to unrealistically 
high values downstream of the reattachment point where a local 
equilibrium state was expected to be almost restored. The reason for 
this behavior is believed to be an underestimation of the turbulent 
viscosity length scales in the inner part of the boundary layer as 
pointed out by Johnson (1985). Since o(x) is the measure of 
nonequilibrium effects, it was initially expected to be close to unity 
after the reattachment point. However, the use of o(x) = 1.0 after the 
reattachment point resulted in the decay of the maximum Reynolds 
stresses being too rapid after the reattachment point. Therefore, 
Johnson suggested that the value of a(x) be limited by an upper bound 
such as 3.0. However, such a limit was not used in the present 
calculations because the value of o(x) did not grow as fast as in the 
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Navier-Stokes solutions of Johnson. It never exceeded a value of about 
3.5. 

The velocity profiles predicted by the Johnson-King and Cebeci- 
Smith models for M = 0.875 are compared with the measurements and 

o© 

Navier-Stokes solutions obtained by the Johnson-King model in Figures 40 
and 41. The agreement between predictions of the Johnson-King model and 
the measurements is especially good in the separated flow region. After 
the reattachment point, the predictions of the Cebeci-Smith model agree 
better with the measurements. This slow flow recovery downstream of the 
reattachment point predicted by the Johnson-King model appears to be a 
shortcoming that requires improvement. However, the effect of this 
discrepancy on the prediction of the surface pressure distribution is 
small. The difference in the predicted velocity profiles obtained by 
the FLARE approximation and windward differencing scheme is negligibly 
small. Figures 42 and 43 compare the profiles of Reynolds shear stress 
calculated by the present method with the measurements and the Navier- 
Stokes solutions obtained with the Johnson-King model for M = 0.875. 

oo 

Overall, the predictions of the Johnson-King model show better agreement 
with the measurements than predictions of the Cebeci-Smith model, 
especially downstream of the bump-sting juncture. However, in the outer 
region of the boundary layer, the Reynolds stresses predicted by the 
Johnson-King model decrease too slowly. 

Figure 44 presents a Mach contour plot obtained by the present 
simultaneous interaction method with the Johnson-King model at 

= 0.875. Note that the plot does not take into account viscous 
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FIGURE 41. Mean velocity profiles for a transonic turbulent flow over 
bump (2) (M = 0.875) 




FIGURE 42. Reynolds shear stress profiles for a transonic turbulent 
flow over a bump (1) (M = 0.875) 








regions explicitly because it was made based only on the inviscid part 
of the interacted solution. The body surface shown in this figure does 
not represents the actual body surface but the shear- layer coordinate 
used in the present calculation. The predicted shock location is in 
good agreement with the measurements as shown in the surface pressure 
comparison. The experimental observations and the Navier-Stokes 
solutions by Johnson (1985) indicated that the shock is highly curved in 

the inviscid region so that the shock wave far from the body surface is 

nearly aligned with the bump-sting juncture. However, the present 
solution shows almost vertical shock wave formation. This is believed 
to be due to the error caused by the velocity potential formulation for 
the inviscid flow. The source of the error is likely the neglect of the 

entropy rise across the shock in the potential flow analysis. To be 

more accurate, the inviscid analysis should satisfy the Rankine-Hugoniot 
relation rather than the isentropic relation across the shock wave and 
the effect of the rotational flow generated by the curved shock should 
be accounted for. 

In the early stage of the present study, the effect of changing the 
modeling constants of the Johnson-King model, a^ and was evaluated 

using the simultaneous interaction method. As discussed by Johnson and 
King (1985), the interaction solution was found to be insensitive to the 
choice of On the other hand, the effect of varying a^ upon the 

prediction is noticeable. Some of the present calculations were 
performed using several different values of a^ The typical value 
deduced from experimental measurements for compressible flows ranges 
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between 0.2 and 0.3. Figures 45-48 compare predictions with different 
values of a. in the range 0.25 ~ 0.4 for the flow at M = 0.875. Note 
that all the Johnson-King model results shown in the previous figures 
were obtained with = 0.25. 

As can be seen in these figures, the results upstream of the shock 
wave seem insensitive to the value of a^, but differences become evident 
downstream of the shock wave. Figure 45 shows that the shock position 
moves very slowly in the downstream direction and the pressure recovery 
in the bump trailing edge region increases with increasing values of a^. 
Interestingly, predictions of the surface pressure with the larger a^ 
(0.35 ~ 0.40) are in better agreement with the measurements of the 2x2 
foot wind tunnel, while the prediction with a smaller value (0.25) of a^ 
agrees much better with the measurements of the 6x6 foot wind tunnel. 

The predicted location of the separation point does not seem to be 
very sensitive to the value used for a^. The location of the 
reattachment point, however, does exhibit a dependence upon a^. With a 
larger value of a , the predicted reattachment point moves upstream. As 
shown by Figures 47 and 48, increasing the value of a^ from 0.25 to 0.35 
appears to give improved predictions for the displacement thickness and 
the maximum Reynolds stresses distributions. However, considering the 
fact that these measurements were taken in the 2x2 foot wind tunnel 
and recalling the observation made in the comparison of the surface 
pressure, these variations are thought to be within the range of errors 
associated with the experimental data. Therefore, it is difficult to 
conclude which value of a^ gives the best overall predictions. It seems 
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FIGURE 46. Skin-friction coefficient distribution with varying a^ for 
transonic turbulent flow over a bump (M = 0.875) 
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FIGURE 47. Displacement thickness distribution with varying a 1 for 
transonic turbulent flow over a bump (M = 0.875) 
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FIGURE 48. Maximum Reynolds shear stress distribution with varying 
for a transonic turbulent flow over a bump (M — 0.875) 
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clear that values of a^ in the range 0.2 ~ 0.3 give reasonably good 
predictions for most of the turbulent flow calculations in this study. 
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VI. CONCLUSIONS 
A. Concluding Remarks 

Two new viscous -inviscid interaction procedures have been developed 
and evaluated for transonic flows. Both employ an inverse finite- 
difference solution of the boundary- layer equations and a direct finite- 
difference solution to the conservative form of the full potential 
equation. The boundary- layer equations are solved in a fully implicit 
manner using Newton linearization with coupling of the continuity and 
momentum equations. The solution procedures have been developed using 
both the SLOR and AF2 schemes to solve the full potential equation. The 
two interaction procedures employ different coupling algorithms; one is 
similar to the procedure developed by Carter (1979) (the semi-inverse 
scheme) and the other, known as the simultaneous methods, is new as 
! applied to finite-difference forms of the boundary- layer and full 

i 

potential equations for transonic flows. The Cebeci-Smith and Johnson- 

I King turbulence models were used for turbulent flow calculations. The 

r 

present schemes were applied to predict flows over three different 
configurations in which significant flow separations occurred; 
incompressible laminar flow over a flat plate with a trough, transonic 

I 

flow over an axisymmetric circular-arc boattail, and transonic flow over 
an axisymmetric circular-arc bump attached to a circular cylinder. From 
the study described in the previous chapters, the following conclusions 
can be made. 
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1. A modification to the semi-inverse interaction method proposed 
by Carter (1979) has been introduced and resulted in a reduction of 
30 ~ 50% in the number of global iterations required for a specified 
level of convergence for the cases considered in this paper. The semi- 
inverse method is found to be stable for calculations of the subcritical 
flow with separation. However, it becomes less stable in supercritical 
flow with increasing size of the separated region so that a smaller 
value of the relaxation factor is required to obtain convergence. The 
value of the relaxation factor required seemed to depend upon the 
turbulence model used. The Johnson-King model required the use of a 
smaller relaxation factor than was needed for the Cebeci-Smith model and 
a larger number of global iterations was generally required for 
convergence. This apparent sensitivity to turbulence models may 
actually be only a sensitivity to the viscous solution being obtained, 
which was significantly different for each of the two models. The 
Johnson-King model tended to predict a larger separated region which 
could be the major cause for the requirement of a smaller relaxation 
factor and more global iterations for convergence. In attempting the 
calculation of the bump flow at M >0.9, the semi-inverse method failed 

oo 

to provide convergence regardless of the turbulence model used. 

2. A new simultaneous interaction method has been developed by 
obtaining the solutions of the boundary- layer equations and the SLOR 
procedure of the potential equation simultaneously. The predictions 
from this interaction method are in good agreement with those from the 
semi-inverse method when the latter gives a converged solution. For the 
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calculation of subcritical flows with a relatively small separated 
region, the simultaneous interaction method is usually about twice as 
efficient as the semi-inverse method. Additionally, the simultaneous 
interaction method was found to be stable over a wider range of flows 
than the semi-inverse method. Also, the simultaneous interaction method 
seems to be relatively insensitive to the turbulence model used. 
Moreover, by employing a pseudo-time dependent approach in the iterative 
procedure of the Newton linearization method for the boundary- layer 
equation, the simultaneous method takes only 20 ~ 30% of the total 
computing time required by the semi-inverse method for supercritical 
flow calculations. 

3. The predictions for turbulent flows are quite dependent upon 
the turbulence model used. Of the two models evaluated, the algebraic 
Cebeci-Smith model generally predicted very poorly for separated flows 
at both subcritical and supercritical Mach numbers. The present 
solutions with the Cebeci-Smith model significantly overpredicted the 
pressure recovery and underpredicted the size of the separation bubble 
and the displacement thickness. The position of the shock wave is 
predicted reasonably well by the Cebeci-Smith model. The turbulence 
model proposed by Johnson and King (1985) was found to provide generally 
better predictions than the algebraic Cebeci-Smith and more 
sophisticated two-equation models for the transonic separated flows 
considered. The predicted locations of the shock wave, separation and 
reattachment points using the Johnson-King model agree very well with 
the experimental data. The displacement thickness is underpredicted in 
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the largely separated flow region and the flow recovery toward the 
equilibrium stage after the reattachment point is also underpredicted by 
the Johnson-King model. The effect of varying the a^ parameter in the 
Johnson-King model was investigated. The value of = 0.25 (used by 
Johnson and King (1985)) is found to give generally satisfactory 
predictions; however, for each separated case computed, it was possible 
to obtained slightly improved predictions by adjusting the a^ value. 
Clearly, the Johnson-King model provides much better predictions than 
the algebraic Cebeci-Smith model for these flows. Yet some room for 
improvement in the turbulence modeling for separated flows remains. 

4. The predictions obtained by the FLARE approximation showed good 
agreement with those obtained by the windward differencing scheme. 

B. Recommendations for Future Study 

The present viscous-inviscid interaction method generally provides 
good predictions for flows considered in the present study. However, 
improvement appears possible in several areas in order to provide better 
predictions for flows with strong interaction. Although the velocity 
potential formulation is accurate enough for most of the transonic flow 
regime, errors associated with the isentropic assumption for the shock 
jump conditions increase with increasing strength of the shock wave. 

This error can be completely eliminated by using the full Euler 
equations, but then computational efficiency must be sacrificed. 

Instead, the full potential formulation can be modified by using an 
entropy correction approach suggested by Klopfer and Nixon (1983) and 
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Hafez and Lovell (1983, 1984). This approach is expected to be 
effective in reducing the error caused by the isentropic assumptions. 

The present simultaneous interaction method is based on the SLOR 
solution procedure for the inviscid flow region. Improvement in the 
rate of convergence might be possible by employing a faster iterative 
relaxation algorithm like the ADI or AF2 schemes as in the case of the 
semi-inverse interaction method. Recently, significant progress has 
been made in improving the efficiency of solution algorithms for the 
time-averaged Navier-Stokes equations. Work should continue on Navier- 
Stokes algorithms and these, too, should be useful in predicting complex 
transonic flow with strong viscous-inviscid interaction. 

The newly proposed Johnson-King model is very simple to use and 
shows great promise for accurately predicting flows with large separated 
regions. However, improvement appears possible. A better description 
of flow recovery toward the equilibrium state, and better modeling of 
the diffusion terms in the transport equation for the maximum Reynolds 
stress are two areas where further work should be done. 
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IX. APPENDIX A: COEFFICIENTS IN THE FINITE-DIFFERENCE 

REPRESENTATIONS OF THE CONTINUITY AND MOMENTUM EQUATIONS 

The coefficients in the finite-difference representation of the 
continuity equation, given by Equation (2.114), are 

a j = 2 Ati -S 

d. = b. 

J J 

h. = 0 
J 


Sj - - 2 A "- (F j + F j-1> 




When the FLARE approximation is used, the coefficients of the 
finite-difference representation of the momentum equation, given by 
Equation (2.115), are 


j An_(K+i) 'V ‘ An + Vi ] 


A. = 


B j = Ad^ry - *r i) - a^ n m ] 


C ~ ~ 9 ~ 

c j = aT ? j (3F j " 2F j } + 2 § B(2F j ' G j) 


~2 


[ -1 1 


A ^mr ^j - - V + K ^ F j - F j-pj 
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D = ' (A. + B ) + 4 |'-S 2 '2F. - F*" 1 ) + 2s 2 BF. 


E j - ttZfW. &Vi - V + K(? j 


H. - i^CG. - ?/) 


i-1 1 ~ 

. J- \ r 1 /n 


s j ■ iUSmr<*i • >¥ F j+1 - V + K(F j - F j-1 


)] 


9 r ~ ~ ^*2 

+ ib F j (F i - v + 2sP(G j - F J > 


When the windward differencing is used, the coefficients of the 
finite-difference representation of the momentum equation, given by 
Equation (2.115), are 


A j An_(K+i) " *j } An/j+i 1 


B j - i^Wj - v - 


C. = ^f-rF.CZFf 1 - 3F ) + 2g p(2F . Z - G ) 
J J J J J J 


+ & 


,i+lwl,; 


An 


& (2x1) - xl> M ± (F - F ) + K(F . - F ,)1 

_(K+l)A$ + UV j V j n K U j+l j' ^ j j-l ;J 


»j ■ - + v - - f 5 +1 > + ^ 




E j - inw)ii; ¥Vi - F j> +K(F j - 
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2 ~~ i+1 ~ ~ ~ 9 

' - V + 2sB(G j - F j > 

Note that the provisional values with tilde (~) are always evaluated at 
ith streamwise station. 
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X. APPENDIX B: COEFFICIENTS IN THE FINITE -DIFFERENCE REPRESENTATION 

OF THE ENERGY EQUATION 

When the FLARE approximation is used, the coefficients in the 
finite-difference representation of the energy equation, given by 
Equation (2.122), are 


A j An_(K+l) [ KAT (,/ 'j " V " An + N l,j+* ] 


B. = 


At _ (K+l) [K Af3 ( *j ' *j } 


AT N U-i ] 


r = c V p i-1 , 2_ 

J A£_ S F j G j + At (K+D^.j+i " ”2 , j _ 2 J 


:(N 0 - N„ , ,) 


D. = - (A. + B.) + T§-g 2 F* 
J J J A? J 


When the windward differencing is used, the coefficients in the 
finite-difference representation of the energy equation, given by 
Equation (2.122), are 


A _ 1 r g /t i ( i+l. 2 v . , 

A j At _ (K+l) [ KA£ + ( *j ‘ *j } " An + N l,j+i ] 


b j ■ ^ 


C j ’ A^ + S F j G j + At_(K+l) (N 2,j+i " N 2,j-i ) 


D. = - (A. + B.) - JLfr] 
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XI. APPENDIX C: RESULTING COEFFICIENTS BY THE USE OF THE MODIFIED 

THOMAS ALGORITHM 


A system of block tridiagonal equations. Equation (2.129), is 
reduced to a set of bidiagonal equations, Equations (2.131) and (2.132), 
by using the modified Thomas algorithm. The coefficients for 2 < j < NJ 
in Equations (2.131) and (2.132) are given as 

* A j 

j Q 0 


C. = ±-[C. - C. . (B. + E.b.) + E.(c. - c. ,)] 
J Q 2 J j-i J J J J J j-i 1 


* _ 1 


H. = r-[H . - H. (B . + E.b.) + E.(h. - h. ,)] 
J Qo J J-l J J J J J J-l J 


S. = ±-[S. - S. (B . + E.b.) + E.(s. - s .' ,) 
J Q 2 J J-l J J J J J j-i^ 


a. = A .Q 
J J 1 


c. = C.Q, + C. ,b . - c. + c. , 
J J 1 J-l J J J-l 


H.Q, + H. ,b, 
J 1 J-l J 


h . + h . 1 

j j-i 


S.Q. + S. b. 
J 1 J-l J 


5 . + S . 

J J-l 


Q - d. + A . b. + a. 

i j j-i j j-i 


Q 2 " “j + + Vl 


The corresponding values at j = 1 are given by Equation (2.133) 
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XII. APPENDIX D: TRANSFORMATION FORMULAS USED FOR THE BOUNDARY 

CONDITIONS IN THE NUMERICAL GRID GENERATION 


The solution procedure for the numerical grid generation was 
presented in Section III. C. It required a Dirichlet type boundary 
condition at all boundaries of the computational domain. In the 
streamwise direction along the body surface, the grid points are 
controlled by the transformation formula suggested by Roberts (1971). 


,, . slnhl “l (£/t max ~ VI , 
X X c sinhfa^a,,) 


where 


“2 = 2 ^ lnI 


1 + (e - l)x c 
•°1 

1 + (e + 1 ) x 


0 < < <*» 

More points are clustered near x as the stretching parameter a becomes 

c 1 

larger. Values of closer to zero result in the grid distribution 
being close to the equal spacing. The typical value of a was between 2 
and 5. At the outer boundary, the grid spacing was set uniform. 

In the normal direction, grid points were concentrated toward the 
wall by the use of an exponential stretching type transformation given 
as 
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<« 3 + 1 ) - (« 3 


“3+1 1 n 

1)( 7) " SX 

OCo “ 1 


r = 


r + 
o 


1 — 

«o + 1 T1 

,3 max + J 

3 


1 < < °° 


This stretching transformation clusters more points near r = as the 
stretching parameter a 3 approaches unity. A value of much larger 
than 1 produces a mesh with approximately uniform spacing. In this 
study, the typical value ranged from 1.05 to 1.1. 
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XII. APPENDIX E: COEFFICIENTS IN THE FINITE -DIFFERENCE 

REPRESENTATIONS USED IN THE NUMERICAL GRID GENERATION 


The Laplace equation which is used to define the coordinate 
transformation is discretized by using a second-order-accurate central 
differencing scheme. The finite-difference representation of the 
Laplace equation in the ADI scheme constructs two sets of tridiagonal 
matrix equations for x and r, Equations (3.41) and (3.42). Since the 
coefficients for x and r are almost identical, only the coefficients for 
x are described explicitly. 

By rearranging the coefficients of the unknown variables. Equations 
(3.41) and (3.42) can be rewritten as 
step 1 


BX . f 1 ? , . + DX.f n . + AX.f n 


i i-l»j 


i i+l,j 


( 12 . 1 ) 


BY.C n . + DY.C n . + AY.C n 

J i.J-1 J i,J J i.J+1 


( 12 . 2 ) 


where 


~n n+1 n 

C . . = x . . - x . . 

i>J i,J 


The coefficients in Equations (12.1) and (12.2) for 2 < i < NI-1 


and 2 < j < NJ-1 are 


AX. = BX. = -A 1 . 1 . 
1 i i, J 


DX. = a + 2A . 
1 
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CX. = awLx n . 
i i.J 


AY. = BY. = -C n . 
J J i,J 


DY. = a + 2C n . 
J i.J 


CY. = f n . 
J i,J 


At the surrounding boundaries of the computational domain, the 


Dirichlet boundary condition is prescribed. Therefore, the coefficients 


in Equations (12.1) and (12.2) at the boundaries are 


AX, = BX, = CX = AX ktt = BX ktt = CX KTT = 0 
1 1 1 NI NI NI 

DX, = DX ktt = 1.0 
1 Ni 

AY 1 * BY 1 = CY 1 = AY NJ - BY NJ “ CY NJ = 0 
DY 1 = DY NJ = 1 ° 


In the above expressions, A£ and An are omitted because A£ = An “ 1.0 
everywhere. The coefficients in the matrix equation for r is easily 
obtained by simply replacing the f^ and x with g^ ^ and r, 
respectively, from the expressions for the above coefficients. 
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XIII. APPENDIX F: COEFFICIENTS IN THE FINITE -DIFFERENCE REPRESENTATION 

OF THE POTENTIAL EQUATION OBTAINED BY THE SLOR SCHEME 


The resulting coefficients in the tridiagonal matrix equations, 
Equation (3.81), obtained from the full potential equation by using the 
SLOR scheme are given as follows. 

For 2 < j < NJ-1 at ith streamwise station, 


A. = - [A2. .5 f (pr) . + R.] 

j i,j r hi j J 


B. = A 2 . . 6 r (pr) - R. . 

j i,j K 'i+$ j-1 


D. = R. + R. + R. + R. , 
J i i-l J J-1 


E. = b.C n . . . + d.C n , . + a,C n , ... + uL 0 n , 

J J i-l, J-1 J i-l, J J i-l, j+1 i,J 


a. = - WCpr,.^ .. ♦ (Pt) i>1+1 i2 lij+1 l 

b j = " [c?r) i-4,j A2 i-i,j + 

d j ' “t R i-r A 2 i.A ( f r >i >J+ 4 1 


where 


1 ^2 

A2 . . = 7 . . 

i,j 4 J 


(13.1) 


and and are given by Equation (3.79). 

At the outer boundary, j = NJ, the potential is fixed by the 


freestream value. Therefore, ^ j = 0 . The corresponding 


coefficients are given as 
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S NJ B NJ E NJ 0 


d nj = 1 -° 


At the wall surface, j = 1, the transpiration velocity boundary 
condition is prescribed. The resulting coefficients are 


A = - 0 R. 
1 r 1 


B 1 = ° 


D, = A4 . + A4 . + 0 R 

1 i l-l r 1 


E. = d C n + a.C n + wL 4? 

1 1 i-l,l 1 i-l,2 i,l 


a l = ' W0 r ( P r) i,l+i A2 l,2 

d l = W t A4 i-l ' 0 r ( P r) i,l+i A2 i,l ] 


A4 = [(pr)4± - ^j)) 

1 J V i+i,l 


0 = l + — 

r r i,l+i 


A closer look at L <t> . reveals the new relation between E and 

1,1 1 

dm /d£ as follows. 


where 


•p i / dm -v . / dm ^ , /■ dm . 

E 1 “ e 0 e l ( d£ h-i e 2 ( d| -) i e 3 ( d£ h+± 


(13.4) 
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where 


e 0 =d l C i-l,l +a l C i-l,2 


+u[5,*4.?. + 9 (pr). ,.,(1 + E +1 )6, + 9 R,f 1 p" , 

Z 1 Z r r 1,1+i n z r 1 7 i J *1,1 


e = - wA5 . .( £ ). i . 
1 l-l p'i-i,l 


- - 2u , 

r i,l ” 1 > 1 


e = wA 5.( £ ) ... . 
3 l p i +f , 1 


and 


A 5 i = ( r } 

A 3 i+i,l 


(13.5) 


By eliminating the upper diagonal terms using a standard Thomas 
algorithm, Equation (3.81) becomes the bidiagonal recurrence 
relationships, given by Equation (3.82). The coefficients of Equation 
(3.82) are given as 


JL. . 1 . 


p. = E. / D. 
J J J 


q. = - B. / D. 
J J J 


D. = D. - A . B . / D 
J J J J+l J+l 


E. = E. - A .E / D 
J J J J+l J+l 
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As given by Equation (3.83), the coefficient p 1 can be more specific in 

.t. 

terms of dm /d£ using the relation given by Equation (13.4). The 
coefficients of Equation (3.83) are then given as 




p 10 l e 0 " A 1 E 2 ^ D 2^ ^ E 1 


P 11 = e l 1 D 1 P 12 = e 2 ' D 1 p 13 = e 3 / D i 
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XIV. APPENDIX G: COEFFICIENTS IN THE FINITE -DIFFERENCE REPRESENTATION 

OF THE POTENTIAL EQUATION OBTAINED BY THE AF2 SCHEME 


The AF2 scheme consists of a two-step procedure. When the main 
flow direction is along the positive £ direction, coefficients in 
Equations (3.88) and (3.89) are given as 


b . = -R . 

J J-l 


d. = (a + R.) 
J J 


e . = awL0 . . 
J i.J 


A. = -R. 

l l 


B. = - (aX + R. ,) 

l 1-1' 


D. = a + aX + R . + R. , 
i l l-l 


E . = f n . + aC n ... 
i i.J i>J+l 


When the flow is in the negative £ direction, A^ and change based on 
the windward differencing scheme. In such a case, 


DftTt 



